Step size is effectively zero

Question

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!

Answer 1

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
*)