# Solution of a BVP by shooting method

I have the next system of ODEs: $$\frac{1}{2}\sigma^2 p''_n(x)+x p'_n(x)+p_n(x)-\frac{p_n(x)-p_{n-1}(x)}{\Delta t}=0,$$ $$p_1(a)=p(b)=0;\quad p_0(x)=0.1,\quad n=1.\dots N,$$ which was derived from the elliptic pde by method of lines: $$\frac{\partial p}{\partial t}=p+x\frac{\partial p}{\partial x}+\frac{1}{2}\sigma^2\frac{\partial^2 p}{\partial x^2},$$ $$lim_{x\rightarrow\pm \infty} p(x,t)=0,\quad \forall t.$$ Now I should set the border values: $a$ and $b.$ I want them to be as far from zero as possible. But the problem arises when I try to solve the system by method of shooting, that doesn't permit me to set them different from $a=0$, $b=1$.

n = 5;
h = 1/n;
U[x_] = Table[Subscript[u, i][x], {i, 0, n}] /.
Subscript[u, 0][x] -> g[x];
Eq = Table[Subscript[eq, i], {i, 0, n}];
\[Sigma] = 0.1;
a[x_, t_] = 1/2 \[Sigma]^2;
b[x_, t_] = x;
c[x_, t_] = -1;
d[x_, t_] = 1;
f[x_, t_] = 0;
g[x_] = 0.1;
\[Alpha][t_] = 0;
\[Beta][t_] = 0;
xl = -0.5;
xr = 0.5;
A=0;
B=1;
eqs = Table[Eq[[i]] =  1/(xr - xl)^2 a[x, h i ] D[U[x][[i]], {x, 2}] + 1/(xr - xl) b[x, h i ] D[U[x][[i]], x] - c[x, h i ] U[x][[i]] - d[x, h i ]  ((U[x][[i]] - U[x][[i - 1]])/h) - f[x, h i ] == 0, {i, 2, n + 1}];
bcs = Table[{U[A][[i]] == 0, U[B][[i]] == 0}, {i, 2, n + 1}]
sols = First[ NDSolve[{eqs, bcs}, U[x], x, Method -> {"Shooting",
"StartingInitialConditions" -> Table[(D[U[x][[i]], x] /. x -> 0) == 1, {i, 2, n + 1}]}]]


When I set $A=-1$ and change initial conditions for method of shooting I get the error:

NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly.

NDSolve::berr: The scaled boundary value residual error of 5.277001462096112*^45 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.

I've also tried to rescale the system by inputing scale variables $xl$ and $xr$. But the same problem holds true. Only current values can be applied. Have no clue what's the problem can be - the analytical solution for the given pde exists uniquely.

• That's just a warning. Did the solution "match the boundary conditions poorly" or not? – Michael E2 Apr 7 '18 at 12:16
• No. The solution bursts at the right border. I've forgotten to attach another message to the question. Now, it's done. – Artem Zefirov Apr 7 '18 at 13:18
• So you are trying to use Method of Lines to discretize in $t$? Usually you discretize in $x$ instead, leaving a set of first order ODEs. – KraZug Apr 9 '18 at 5:30
• Some authors assert that discretization in t is possible and is better in some extent. Cause we are capable to make use an analytical solution for BVP having infinitive domain. In the book written by Gunter H. Meyer the discretization in t is used to deal with the diffusion equation in a finite domain. – Artem Zefirov Apr 9 '18 at 9:04
• I don't know much about numerical methods, just what I've done before. – KraZug Apr 9 '18 at 10:32

The main issue with your code is very simple, you are trying to take your domain in $A\leq x\leq B$ and split it into $n$ equal sized pieces. These should have size $h=(B-A)/n$, but you are setting them to be size $h=1/n$. All you need to do for your code to not give errors is to define A and B first, and then set h=(B-A)/n.

It does seem very slow though for larger values of n, and you might want to use higher order derivative formulae.

Adapting some code I have for method of lines, but discretizing in $x$ instead of $t$ (have to take a lot of points here with the default settings, else error creeps in and the solutions blow up.):

A=-5;B=5;σ=0.1;
n=300;h = (B-A)/n;
gridpts=(Range[A,B,h]);
P[t_] = Table[Subscript[p,i][t],{i,0,n}];
nDeriv[m_?NumericQ,fn_]:=NDSolveFiniteDifferenceDerivative[m,gridpts, fn,DifferenceOrder->2];
psubs={p[x,t]:>nDeriv[0, P[t]],Derivative[Pattern[i,Blank[]],0][p][x,t]:>nDeriv[i,P[t]],Derivative[0,Pattern[i,Blank[]]][p][x,t]:>D[P[t],{t,i}],XX[x]:>gridpts};
eqn=D[p[x,t],t]==p[x,t]+ XX[x] D[p[x,t],x]+1/2 σ^2 D[p[x,t],x,x];

(* find an initial perturbation from zero satisfies the boundary conditions *)

initsol=(a0+a1 x+a2 x^2)/.Solve[{p[A]==0,p[B]==0}/.p->Function[{x},(a0+a1 x +a2 x^2)]][[1]];

lines=NDSolve[{MoLeqns,inits},Table[Subscript[p,i],{i,0,n}],{t,0,50}][[1]];
tend=lines[[1,2,1,1,2]];


This gives a steady solution that has a spike at $x=0$ and rapidly decays outside that:

ListLinePlot[P[tend] /. lines, PlotRange -> All,  DataRange -> {A, B}, AxesLabel -> {"x", "p"}]


Or you can see how the individual lines behave:

ParametricPlot3D[Evaluate[Table[{A+ i h,t, Subscript[p,i][t] /.lines},{i,0,n}]],{t,0,tend}, PlotRange->All, AxesLabel->{"x","t","p"}]


• I've fixed the error concerning h but it didn't help. I've tried also to use the second order derivative in t, but it doesn't fix the problem. Moreover, application of the second order derivative with current settings produces the next warning: NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. If a solution is computed, it may match the boundary conditions poorly. – Artem Zefirov Apr 9 '18 at 9:17
• Hmm, If I set h as defined then I can move B away from 1, e.g. 20. But trying to reduce A doesn't work. – KraZug Apr 9 '18 at 10:26
• Thank you. But it still interesting what's wrong with discritization in t. Besides, have you any idea why does the stretching of the domain lead to incorrect results. I've set A=-20, B=20. Can normalization fix it? – Artem Zefirov Apr 9 '18 at 15:58
• For my code, with your particular equations you'll need a finer discretization (larger n) or better error control (options for NDSolve) if you try and increase the size of the domain. You can see that oscillations start to grow if you are not careful. But looking at the plot it is very localised, so $\infty =5$ should be fine. – KraZug Apr 9 '18 at 19:07