# For loop to perform a shooting search in NDSolve

First comes the disclosure: I am very new to Mathematica and not so much a numerical guy in general. I am having a hard time using the for loop routine in Mathematica the correct way as well as finding resources as how to do so.

I heard that we should avoid as possible using For loops in the first place in Mathematica but I don't see how I could use an alternative to solve my problem and I do not manage to make my for loop working.

Here is my problem: I want to find a way to search for an undefined boundary condition according to some smoothness condition by varying one of the initial condition. In my example I want to find the point for which both solution reaches zero at the same time. The algorithm should be able to

1. Make a guess for X0
2. Compute the Solution
3. If X reaches 0 before Y, increase the guess, otherwise decrease it
4. iterate until X(T) = Y(T) = 0

I manage to do 1 2 3 and but do not manage to get to update on my guess. Here is the code I was writing:

crit = 0.01;
f[x_] := If[Abs[x] < crit, 1, 0];

For[ite = 1; ql = 0; qh = 5; F0 = (ql + qh)/2; a = 100,

ite < 100; f[a] > 0,

ite++; ql = qlu; qh = qhu,

sol = NDSolve[{y'[t] - 2 y[t] == 4 - t, x'[t] == 2 x[t] + 4 - t,
y == -1.5, x == F0,
WhenEvent[{x[t] == 0, y[t] == 0}, end = t;
"StopIntegration"]}, {x, y}, {t, 0, 20}];

xend = x[end] /. sol;
yend = y[end] /. sol;
a = First[xend] - First[yend];
If[Abs[First[xend]] < crit, qlu = F0, qhu = F0];
Print[qlu, qhu]

]


My two questions: 1. How do I make Mathematica update my guess ? 2. Am I right to use a For loop ? What could be an alternative (better) option ?

• I guess it wouldn't hurt to make the For loop work in this case so that you can clarify your logic. It never starts because f=0. If you fix that you'll throw an error because end is never set (unless your initial guess is the solution). A better approach would be to use FindRoot – george2079 Oct 22 '15 at 19:42
• Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Dr. belisarius Oct 22 '15 at 19:46
• WhenEvent[{x[t] == 0, y[t] == 0} .. stops when x OR y is zero. If you want both do WhenEvent[x[t] == 0&&y[t] == 0,.. (I suspect your system has no solution though ) – george2079 Oct 22 '15 at 19:51
• make that no non-trivial solution..:) – george2079 Oct 22 '15 at 20:45

y'[t] - 2 y[t] == 4 - t,
x'[t] == 2 x[t] + 4 - t


... they are equal!

So for getting x and y going to zero at the same time just do x = y ...

Edit

Now,more seriously, I believe you mis-engineered your example. Something more general (although not completely) is:

Please note that the functions are coupled now:

k[x0_?NumericQ] := Module[{res, s, f0},
s = ParametricNDSolveValue[{
y'[t] == 2 y[t] x[t] + 4 - t,
x'[t] == 2 x[t] + 4 - t,
y == -1.5, x == f0,
WhenEvent[{x[t] == 0}, res =  y@t^2; "StopIntegration"],
WhenEvent[{y[t] == 0}, res =  x@t^2; "StopIntegration"]},
{x, y}, {t, 0, 20}, f0];
s[x0];
res];

Plot[k[f0], {f0, -2, 0}] And we find the crossing as:

FindRoot[k[f0], {f0, -2}]
(* -0.907017 *)


Please note that the res values have been defined so that we can use FindRoot[ ] .... in this answer's history record you may find other approaches :)

edit 2

Somewhat more elegant

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
s = ParametricNDSolveValue[
{y'[t] == 2 y[t] x[t] + 4 - t,
x'[t] == 2 x[t] + 4 - t,
y == -1.5, x == f0,
WhenEvent[{x[t] == 0, y[t] == 0}, "StopIntegration"]},
{x, y}, {t, 0, 20}, f0];
tf[t_] := InterpolatingFunctionDomain[Last@s[t]][[-1, -1]]
val[t_?NumericQ] := Subtract @@ (Through[s[t][tf[t]]]^2)
FindRoot[val[f0], {f0, -2}]
Plot[val[x], {x, -2, -10^-10}]

(* {f0 -> -0.907017} *) 