2
$\begingroup$

I'have a 2d-system of differential equations, analitycally solved, depending on a parameter. I know that, by continuity, there exist a minimum value of the parameter such that trajectory passes through a specific point. My problem is to find such value (also approximatively).

This is the system-equation:

ti = 0;
yi = 0;
zi = -.75;
zf = -.5;
eps=.01;

sol = FullSimplify[DSolve[{y'[t] == 1/(2 zf) y[t] - u z[t], 
z'[t] == -1 + 1/zf z[t] + u y[t], y[ti] == yi, z[ti] == zi},{y[t],z[t]},t]];

The parameter is u.

I have to find the first u such that trajectorie intersect the point: z=zf; y=Sqrt[-2*(zf)*eps + eps^2]

Note that this point is done by the intersection between z=zf and the circle of radius |zf|+eps. In the following, you can see an animation of the system:

Manipulate[
 Module[{sol, y, z, t}, 
  sol = First@DSolve[{y'[t] == 1/(2 zf) y[t] - u z[t], 
  z'[t] == -1 + 1/zf z[t] + u y[t], y[ti] == yi, 
  z[ti] == zi}, {y[t], z[t]}, t];
Show[p1, p2, Graphics[{Red, Line[{{0, -1}, {0, 1}}]}],
Graphics[{DotDashed, Red, Thickness[.006], 
 Line[{{-1, zf}, {1, zf}}]}],
ParametricPlot[{y[t] /. sol, z[t] /. sol}, {t, 0, tend}, 
PlotStyle -> Thickness[.004]]]],
      {{tend, .1, "t"}, .01, 20, .1, Appearance -> "Labeled"},
      {{u, 10, "u"}, -300, 300, .5, Appearance -> "Labeled"},
      {{zi, -.75, "zi"}, -1, 1, Appearance -> "Labeled"},
      {{yi, 0, "yi"}, -1, 1, Appearance -> "Labeled"},
Initialization :>
(deltaA[y_, z_] := 1/(2 zf) y^2 - z + 1/zf z^2;
circ[y_, z_] := y^2 + z^2;
p2 = ContourPlot[{circ[y, z] == (-zf + eps)^2}, {y, -1, 
  1}, {z, -1, 1}, ContourStyle -> Yellow, GridLines -> Automatic, 
 Frame -> True, FrameLabel -> {"y", "z"}, RotateLabel -> False, 
 LabelStyle -> {FontSize -> 20}];
p1 = ContourPlot[{deltaA[y, z] == 0}, {y, -1, 1}, {z, -1,
   1}, ContourStyle -> Green, GridLines -> Automatic, 
 Frame -> True, FrameLabel -> {"y", "z"}, RotateLabel -> False, 
 LabelStyle -> {FontSize -> 20}];)]
$\endgroup$
  • $\begingroup$ y depend on t, the time, AND u, a free parameter ... $\endgroup$ – Mike84 Jul 16 '14 at 17:44
  • $\begingroup$ I'have just added a Manipulate to show the behaviour of the system. Anyway, maybe I've not understand your question...the definitions of y[t],z[t] come from DSolve. They are analytic.. $\endgroup$ – Mike84 Jul 16 '14 at 18:18
  • $\begingroup$ Maybe I am the one not understanding :) Your code doesn't compute, brackets mismatch. $\endgroup$ – Öskå Jul 16 '14 at 18:21
  • $\begingroup$ ??? :) Try again...If your problem was in DSolve, I'have rewrited it $\endgroup$ – Mike84 Jul 16 '14 at 18:29
  • $\begingroup$ I'm afraid that it doesn't work. $\endgroup$ – Öskå Jul 16 '14 at 18:38
2
$\begingroup$

Here is a stab at what I think you are asking:

 ti = 0;
 yi = 0;
 zi = -.75;
 zf = -.5;
 eps = .01;
 sol = {y[t], z[t]} /. First@FullSimplify[DSolve[{
   y'[t] == 1/(2 zf) y[t] - u z[t],
   z'[t] == -1 + 1/zf z[t] + u y[t],
   y[ti] == yi, z[ti] == zi},
        {y[t], z[t]}, t]];
 target = {Sqrt[-2*(zf)*eps + eps^2], zf};

brute force discretize the solution and find the minimum distance to the target point (for such a highly nonlinear function this is much faster than NMinimise, and guarantees finding a global min, within the discretization approximation of course )

 dis[u0_] := 
     Min@Table[Norm[ (sol /. u -> u0) - target ] , {t, 0, 10, .05}];

note the range and increment on t here are important tuning parameters to play with.

 Plot[dis[u], {u, -1, 1}]

enter image description here

you come very close to your target point around 0.21.

Armed with a good guess now we can use FindMinimum:

 FindMinimum[Norm[sol - target], {u, .21 } , {t, 1}]

{3.23748*10^-9, {u -> 0.230625, t -> 1.62219}}

 ParametricPlot[
    Table[Chop[sol /. u -> u0 ], {u0, {.1, .230625, .3, 1, 5}}] , {t, 0, 10},
        Epilog -> Point[target], PlotRange -> All, AspectRatio -> 1]

enter image description here

$\endgroup$
  • $\begingroup$ thank you very much for this very useful and didactic answer. I have a question: why do you need chop? Furthermore I try to make a numerical procedure, using a Do and the function FindAllCrossing. But it doesn't work and I don't know why. Can I post an answer with my wrong-code to understand my error? $\endgroup$ – Mike84 Jul 17 '14 at 0:04
  • $\begingroup$ Chop removes a numerically negligible complex part. Hitting a point in 2d isn't formulated in terms of 'zero crossings' so i cant see how you could use the methods in the link. You need to find minimum distance to the point and if you need to automate it you need to put in a check that the minimum you find is really zero. $\endgroup$ – george2079 Jul 17 '14 at 12:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.