# Solving a bad-behaving equation with shooting method

I would like to numerically solve

$$a^2\frac{d}{d\sigma}[(1+A_tQ)\frac{d\psi}{d\sigma}]=\psi(1+A_tQ-a\Phi\frac{dQ}{d\sigma})$$

where $$Q(\sigma)=\frac{2}{\sqrt{\pi}}\int_0^\sigma exp(-x^2)dx$$.

$$a$$ and $$A_t$$ are certain parameters and $$\psi(x)$$is the desired function with certain boundary values and $$\Phi$$ is the eigenvalue of the equation. I started with the built-in shooting method:

a = 1/5;
A = 0.843;
Q[x_] := 2/Sqrt[Pi]*Integrate[Exp[-p^2], {p, 0, x}];
eq := a^2*D[(1 + A*Q[x])*D[psi[x], x], x] ==  psi[x]*(1 + A*Q[x] - a*Phi[x]*2/Sqrt[Pi]*Exp[-x^2]);
bc = {psi'[-2] == -1/a*Exp[-2/a], psi[-2] == -Exp[-2/a], psi[100] == 0};
s = NDSolve[{eq, bc, Phi'[x] == 0}, {psi[x], Phi[-2]}, {x, -2, 100},
Method -> {"Shooting",
"StartingInitialConditions" -> {Phi[-2] == 3}}];


However, this equation behaves badly, and with shooting method I always get NDSolve:berr, and the solution of Phi stays where the initial condition is set.

So I wrote a simple shooting method myself:

Q[x_] := 2/Sqrt[Pi]*Integrate[Exp[-p^2], {p, 0, x}];
eq := a^2*D[(1 + A*Q[x])*D[psi[x], x], x] == psi[x]*(1 + A*Q[x] - a*Phi*2/Sqrt[Pi]*Exp[-x^2]);
bc := {psi'[-2] == -1/a*Exp[-2/a], psi[-2] == -Exp[-2/a]};
a = 1/5;
A = 0.843;
Phi = 3;
d = 0.1;
Sol = NDSolve[{eq, bc}, psi, {x, -2, 100},
Method -> "ExplicitRungeKutta"];
m1 = Maximize[{Abs[psi[x] /. Sol[[1]]], 80 <= x <= 100}, x][[1]];
While[d > 0.0001,
Phi = Phi + d;
Sol = NDSolve[{eq, bc}, psi, {x, -2, 100},
Method -> "ExplicitRungeKutta"];
m1new = Maximize[{Abs[psi[x] /. Sol[[1]]], 80 <= x <= 100}, x][[1]];
If[m1new > m1, Phi = Phi - d; d = d/2, m1 = m1new];
Print[Phi];
]


This just works better and gave an approximate value, although the solution of $$\psi(\sigma)$$ is still diverging.

Since my scanning scheme is rough and doesn't guarantee a best solution in principle, is there any way to get the built-in shooting method to work? Or is there any better approach to this job?

• "...and the solution of Phi stays where the initial condition is set" doesn't surprise because of ode Phi'[x]==0! Commented Apr 17 at 10:38
• @UlrichNeumann By this I mean that, the shooting method always returned a value of Phi where my "StartingInitialConditions" is set, 3 in my sample code. However I expect shooting method to return a different optimized value of Phi. I guess this is not related to Phi'[x]==0, as this just constrain it to a constant. Commented Apr 17 at 11:34

## 1 Answer

To solve this problem we don't need shouting method, also range $$-2\le x\le 100$$ is to large since even for $$-2\le x\le 5$$ we should compute with WorkingPrecision -> 30. Numerical solution is given by

a = 1/5;
A = 843/1000;
Q[x_] := Erf[x];
eq = a^2*D[(1 + A*Q[x])*D[psi[x], x], x] ==
psi[x]*(1 + A*Q[x] - a*Phi*2/Sqrt[Pi]*Exp[-x^2]);
bc = {psi'[-2] == -1/a*Exp[-2/a], psi[-2] == -Exp[-2/a]};
s = ParametricNDSolveValue[{eq, bc}, psi[5], {x, -2, 5}, {Phi},
WorkingPrecision -> 30];

sol = FindRoot[s[z] == 0, {z, 378/100}]

(*Out[]= {z -> 3.78066}*)

Phi0 = Rationalize[z /. sol, 10^-31];

sol1 = NDSolveValue[{eq /. Phi -> Phi0, bc}, psi, {x, -2, 5},
WorkingPrecision -> 30];


Visualization

Plot[sol1[x], {x, -2, 5}, Frame -> True, GridLines -> Automatic]


As we can see from this picture we can't extend numerical solution up to x=100 since asymptotic solution psi is proportional to $$e^{-5x}$$, therefore psi[100] is about Exp[-500].

• Thank you for the answer. I chose shooting method because of the original paper I read :Phys. Fluids 5, 417–425 (1962), and it seems 'FindRoot' in your code does the job automatically. But I'm still confused, as the solution is still divergent for large x, however as you mentioned 'psi' should converge to 0 as $e^{-5x}$. In fact limiting the region to '[-2,5]' leads to some different values of 'Phi', if I change the initial trial value in 'FindRoot'. How should I judge from the different results? I would hope to generalize this problem for different values of 'a&A' so this seems to matter. Commented Apr 18 at 5:50
• @Repentanze There is unstable mode as well proportional to $e^{5 x}$ or $e^{x/a}$ in general case. So if we need to describe instability we should combine modes. Commented Apr 18 at 11:27