2

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}];)]
george2079
  • 38,913
  • 1
  • 43
  • 110
Mike84
  • 77
  • 6

1 Answers1

2

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

george2079
  • 38,913
  • 1
  • 43
  • 110
  • 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? – Mike84 Jul 17 '14 at 00:04
  • 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. – george2079 Jul 17 '14 at 12:49