# Step size is effectively zero

I've been trying to solve Bodewadt flow equations which is a system of differential equations.

\begin{align*} F^2 - G^2 + HF' - F'' + 1 &= 0 \\ 2GF + HG'-G'' &= 0 \\ 2F + H' &= 0 \end{align*}

with boundary conditions:

\begin{align*} F(0)=G(0)=H(0)=0\\ F(\infty)=0,\;G(\infty)=1 \end{align*}

I turned them into a system of first order ordinary differential equations where $$F=x;F'=y;G=z,G'=s$$ and $$H=p$$. Thing is, I don't have the derivatives $$F'(0)$$ and $$G'(0)$$, so to make use of the infinite boundary conditions instead, I use a sufficiently large number instead of $$\infty$$. Approximately, $$14$$ works fine.

sol = {x[t], y[t], z[t], s[t], p[t]} /.
NDSolve[{x'[t] == y[t], y'[t] == x[t]^2 - z[t]^2 + s[t] y[t] + 1,
s'[t] == 2 z[t] x[t] + p[t] s[t], z'[t] == s[t], p'[t] == -2 x[t],
x[0] == 0, z[0] == 0, p[0] == 0, x[14] == 0, z[14] == 1}, {x, y,
z, s, p}, {t, 0, 14}]


The problem is that when I run the code I have the next problem:

NDSolve::ndsz: At t == 3.4508216573870163, step size is effectively zero; singularity or stiff system suspected.

I really am not acquainted with Mathematica, so I really don't how to solve a system of ODEs with boundary conditions at infinity.

Thanks!

• there is a tag boundary-conditions-at-infinity for questions related to your problem. It’s almost certainly the ”StartingInitialConditions” that are chosen for the ”ShootingMethod” used for BVPs, which you can also look up in the docs. Finally Method -> “FiniteElement” might work for you. Commented Dec 24, 2020 at 15:23
• Hi, I really do not understand what you mean, I am not a native speaker, I guess BVPs stands for boundary value problems, is it a book? I don't know what you mean by docs, I am not acquainted with this page neither. Thanks! Commented Dec 24, 2020 at 16:54
• "Docs" is short for "Documentation" or "Help" pages, like this one: reference.wolfram.com/language/tutorial/NDSolveBVP.html, which discusses Boundary Value Problems (BVP). You can also get to it through the menu Help > Wolfram Documentation in the desktop version (search for "BVP"). -- If click on the tag in my first comment, you will get a list of posts on this site dealing with boundary conditions at infinity. Commented Dec 24, 2020 at 17:31
• Thanks for explanation! Yet, when i click on the tag boundary-conditions-at-infinity of you first comment it says there no questioned tagged. Anyway, the documentation seems to be helpful. Commented Dec 24, 2020 at 17:46
• Sorry about that. I was on a phone and mistyped. Try boundary-condition-at-infinity Commented Dec 24, 2020 at 17:50

I think that ODE systems with $$\| X'\| \sim O(\| X \|^2)$$ tend to be unstable, that is, a small rounding error has a chance to cause a solution to blow up. First, boundary-value problems (BVPs) are not guaranteed to have solutions, and without a proof or evidence that a solution exists, difficulty in solving one should raise the question whether there is a solution to find. Second, instability makes a BVP at infinity hard to solve numerically, even when solutions do exist.

I haven't made a thorough analysis of the problem, but I got lucky (I think it was luck). When you get a NDSolve::ndsz error and no solution returned when solving a BVP, it is usually because of bad starting initial conditions that yield a singularity, which conditions are automatically chosen for the shooting method. The shooting method computes the solution to an initial-value problem (IVP) and adjusts the initial conditions (ICs) trying to get the solution to converge to the boundary conditions (BCs).

You can control the "Shooting" method of NDSolve somewhat (see https://reference.wolfram.com/language/tutorial/NDSolveBVP.html). You can specify "StartingInitialConditions", but sometimes this isn't enough. In such cases one might use ParametricNDSolveValue to manually implement a shooting method (How to avoid NDSolve::ndsz problem (singularity problem), Nonlinear differential equation: numerical solution).

In this case I tried "StartingInitialConditions". The first guess gave NDSolve::nderr error test failure. If you know your ODE system, you might have insight into what are good guesses and what are not. I didn't have such insight, but I got lucky. To give me a better chance, I reduce the precision and accuracy goals and raised the WorkingPrecision. The goal was to get a pretty good solution and then to refine it.

odes = {x'[t] == y[t], y'[t] == x[t]^2 - z[t]^2 + s[t] y[t] + 1,
s'[t] == 2 z[t] x[t] + p[t] s[t], z'[t] == s[t], p'[t] == -2 x[t]};
bcs = {x[0] == 0, z[0] == 0, p[0] == 0, x[14] == 0, z[14] == 1};
vars = {x, y, z, s, p};
sol32 =
NDSolve[{odes, bcs}, vars, {t, 0, 14},
Method -> {"Shooting",
"StartingInitialConditions" -> {s[14] == 0, y[14] == 0,
p[14] == 0, x[14] == 0, z[14] == 1}}, PrecisionGoal -> 6,
AccuracyGoal -> 6, WorkingPrecision -> 32];

solMP =
NDSolve[{odes, bcs}, vars, {t, 0, 14},
Method -> {"Shooting", "StartingInitialConditions" -> {
Through[vars[14]] == Through[vars[14] /. First[sol32]]}
},
PrecisionGoal -> 10, AccuracyGoal -> 10
]

ListLinePlot[vars /. First[solMP], PlotLegends -> vars]


Check the BCs:

bcs /. Equal -> Subtract /. First[sol32] // Norm
bcs /. Equal -> Subtract /. First[solMP] // Norm
(*
0.0337144277985927458402168128813
5.99878*10^-6
*)

• Thanks! I appreciate your answer. It is, indeed, the solution I had been seeking. Commented Dec 24, 2020 at 18:42
• @SebastiánFrades You’re welcome! Commented Dec 24, 2020 at 19:09
• Luckily, the built-in shooting method worked right out of the box. If it gets bad enough, one would usually have to use ParametricNDSolve(Value)` to manually implement shooting. Commented Dec 25, 2020 at 3:58
• @J.M.'sdiscontentment Yep, it's mentioned in the last sentence of the third paragraph. Commented Dec 25, 2020 at 4:06