Re: NDSolve with a constraint : how ?

*To*: mathgroup at smc.vnet.net*Subject*: [mg72597] [mg72597] Re: NDSolve with a constraint : how ?*From*: dh <dh at metrohm.ch>*Date*: Thu, 11 Jan 2007 03:44:46 -0500 (EST)*References*: <en82l6$6pp$1@smc.vnet.net>

Hi Cham, here is a 2 dim. example for your problem. Imagine that you throw a stone with a angle p and given velocity upward. What is the angle to come as near as possible to a given point (the stone has not enough energy to reach the point)? Method: we define a function that calculates the distance of closest approach and use this function in FindMinimum: v=5; (*velocity*) g=10; (*acceleration*) reference={2,1}; (*point to approach*) f:={fx[#],fy[#]}&; (*general 2 dim function*) dist[p_Real]:=( eq=Sequence@@Thread[#]&/@{ f''[t]=={0,-g}, (*these are Newtons equations *) f[0]=={0,0},(*initial conditions*) f'[0]==v{Sin[p],Cos[p]} }; res=f[t] /. NDSolve[eq,f[t],{t,0,1}][[1]]; (*the trajectory*) t0=FindRoot[ D[Plus@@((res-reference)^2), t]==0,{t,.5}];(*parameter t for closets approach*) di=Plus@@((res/.t0)- reference)^2; (*calculate the distance of the closest approach*) Print["p=",p,", dist= ",di]; ParametricPlot[res,{t,0,1}]; (*make a picture*) di ); FindMinimum[dist[p],{p,0.5}] Daniel Cham wrote: > I need to find a proper way to solve an equation with the NDSolve operation. I'm looking for a solution { x[ t ], y[ t ], z[ t ] } which should obey some constraint (the initial conditions { x[0], y[0], z[0] } are not well known and I only have some very approximate values). How should I do this ? > > More specifically, I'm using a simple code like this : > > NDSolve[ > { > x'[t] == Fx[ x[t], y[t], z[t] ], > y'[t] == Fy[ x[t], y[t], z[t] ], > z'[t] == Fz[ x[t], y[t], z[t] ], > > x[0] == x0, > y[0] == y0, > z[0] == z0 > }, > > {x, y, z}, {t, 0, 100} > ] > > Mathematica then finds easily a solution. But the solution I'm looking for must obey a constraint, and the inital conditions {x0, y0, z0} aren't well known. I need to find the initial values {x0, y0, z0} for which the horizontal distance is the closest to some constant, for an unknown "t" : > > rho = Sqrt[x[t]^2 + y[t]^2] = cste, for some unknown "t". > > How can I program Mathematica so it could find the right set of numbers {x0, y0, z0} ? > > For the moment, all I can do is find a solution from some approximate values {x0, y0, z0}, then check by trial and errors if there's a "t" which gives rho = ctse (approximately). If not, I have to use the NDSolve again, and again (changing a bit the inital values after each trial), until it works. This can be very long, especially since I don't know what value of "t" will satisfy the constraint. > > Any better idea ? >