# invert system of two parametric ODEs and plot results of do-loop

I can solve this system of parametric ODE (u is the parameter,while g,e and zi are constants depends on the physic of the problem)

zf = -1/2;
g = 1;
e = .2;
Y[t_] :=  g/(2 zf) y[t] + u z[t];
Z[t_] := -g + g/zf z[t] - u y[t];
sol =
FullSimplify[
DSolve[{y'[t] == Y[t], z'[t] == Z[t],y == 0 , z == -zf},
{y[t], z[t]}, t]]


Then, I make the assignement:

zsol[t_] := Evaluate[z[t] /. sol[][]]
ysol[t_] := Evaluate[y[t] /. sol[][]]


My problem is to compare the time t(u) taken by trajectories {y[t,u],z[t,u]} to reach a neighborhood of zf whit radius e, with the time trelax, relative to u=0 (zf is the fixed point of free dynamic). Not for all u the problem has solution: in general, for u large enough, trajectories reaches some fixed point in an infinity time following a spiral (as I seen plotting the system for various u), cause of the exponential-complex solutions, and there's a set of u without solutions (i.e, a set of u for which trajectories pass no close to zf). Then I have thought to use an iterative costruct, a do-loop mixed with an If. For simplicity, I've started, in particular, with:

trelax = Log[(-e/(2 zf))]*zf/g;
Do[Do[
If[N[Abs[zsol[T]] <= (zf + e) && N[Abs[ysol[T]] <= e,
If[NSolve[ysol[t], t] <= trelax],
Print[{T, u}]]], {T, 0, .9, .01}], {u, 0, 40, .2}]


Remark_1: trelax is the correct analytic result of dynamical system for u=0. Remark_2: I expect from this code to obtain a list of couple {t,u}, with t<=trelax and u>0, but really I don't obtain nothing! I think I'm wrong in the definition of the variable in the function...and in many other things!

Finally I ask you if is it possible to change the loop to obtain directly a plot {t,u}, i.e. a curve to compare with constant trelax...

Thank you

P.S: If you want to see some plots, I attach the ad hoc code:

Clear[u]
u = 5;
ParametricPlot[{y[t] /. sol[][], z[t] /. sol[][]}, {t, 0, 10}, PlotRange-> All]

• Your Do loops as written have several missing brackets and incorrect capitalisation (Nsolve instead of NSolve). What's the numerical value for zf? NSolve gives a rule as a result, not a value which you can directly compare. Maybe you want Sow and Reap rather than Print. Feb 18, 2014 at 20:21
• ok, thanks! actually I have written zi instead zf in the initial definitions. In my code (on my pc...) NSolve is written correctly, but don't work. Now I try to learn all about Reap and Sow...But, if you have understand my problem, do you think that my idea is right? What do I use instead NSolve? Feb 18, 2014 at 20:46

This was (partially answered on the mathematica-community site).

Here is that posting:

Without solving the problem of finding the condition where the trajectory lies within a disk. This modification of your code may help you construct a numerical technique.

sol = DSolve[{y'[t] == g/(2 zf) y[t] + u z[t],
z'[t] == -g + g/zf z[t] - u y[t], y == 0, z == -zf},
{y[t], z[t]}, t]


This is your original equation, but without putting in the defintions for Z, etc.

Extract the solutions and simplify

{{ysol, zsol}} = Simplify[{y[t], z[t]} /. sol,
Assumptions -> Element[u, Reals] && Element[zf, Reals] && t > 0 ]


Explore the solutions...

With[{y = ysol, z = zsol},
Manipulate[ParametricPlot[{y, z}, {t, 0, tend}],
{{tend, 0.759}, 0, 100},
{{g, 0.354692}, 0, 10},
{{zf, -1.8734}, -5, 10},
{{u, 0.922386}, 0, 10}
]
]


You can use a numerical method to locate the conditions where the trajectory is within a given distance to a point, or use the WhenEvent option with NDSolve

In response to your comment,

I am not sure what you were trying to do, but here is an example of finding a (not necessarily unique) set of parameters that are close to the point (1,1).

FindMinimum[EuclideanDistance[{ysol, zsol}, {1, 1}], {g, u, zf, t}]


If you plug the result into the Manipulate (i.e., type them in by clicking the the "+" next to the slider), you will see that it finds a good approximation to a solution.

• thank you, Craig. You have taken me a pretty suggestion and Manipulate is a very useful tool to have an idea about what is happening. But I'm a novice in mathematica...What do you have in mind when you talk about "numerical method etc..."?: in fact my greatest difficult is how to extract a list from the loop, and plot it... Feb 19, 2014 at 0:30
• @Mike84, I wasn't sure what you were looking for, but here is an example of how to find a trajectory that approaches a neighborhood of (1,1): Feb 19, 2014 at 10:41
• To begin I'm looking for the first u such that, given the other parameters, trajectory passes as close as possible to zf (say, an arbitrary small distance e). This problem may be solved with your methods. But for second step, I wish to find a graph {t*,u*} that compare the time t* taken to reach (a neighborhood of) zf given u*...It's for such problem that I tried with a do-loop (or similar)...Thank you very much! Feb 19, 2014 at 14:53