# How to avoid NDSolve::ndsz problem (singularity problem)

I have the next model

Vs = px 1/YXS (Mm x2[r])/(KS + x2[r]) x1[r]/(KO2 + x1[r])

Vo = px 1/YXO (Mm x2[r])/(KS + x2[r]) x1[r]/(KO2 + x1[r])

eqs = {{x1''[r] + 2/r x1'[r] == Vo/DO2, x1'[R0] == 0,
x1[R] == Cb}, {x2''[r] + 2/r x2'[r] == Vs/Ds, x2'[R0] == 0, x2[R] == Sp}}

R = 0.4; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.2; R0 = 0.000001; YXO = 0.42857;


When I try to use NDSolve solution for any R>0.3 it doesn´t work. Several errors occur:

sol = NDSolve[eqs, {x1, x2}, {r, R0, R}];

NDSolve::ndsz: At r == 0.26951986669987066, step size is effectively zero; singularity or stiff system suspected. >>
NDSolve::ndsz: At r == 0.18524442932777493, step size is effectively zero; singularity or stiff system suspected. >>
NDSolve::ndsz: At r == 0.21451570142533177, step size is effectively zero; singularity or stiff system suspected. >>
General::stop: Further output of NDSolve::ndsz will be suppressed during this calculation. >>
FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option. >>
NDSolve::berr: There are significant errors {0.0000127787,-8.2,5.9963*10^-6,-9.04008} in the boundary value residuals. Returning the best solution found. >>


• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) 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! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Jan 29 '15 at 2:41
• I only get the first error message. Certainly no FindRoot errors -- ??? – Michael E2 Jan 29 '15 at 3:17
• The OP must have V9. I tested it on V10. – Michael E2 Jan 29 '15 at 13:03
• Yeah! As you notice I'm new here. Previously, I found great answers to my problems on this site, without being part of it. Now really had to ask. I will follow your instructions carefully and thanks. – myguelangelo Jan 29 '15 at 17:11
• Yes as you say I have version 9 – myguelangelo Jan 29 '15 at 17:13

Update notice 2: Found starting initial conditions that work in V10.

Update notice: bbgodfrey pointed out that the built-in shooting method does not work in V10.0.1 (nor V10.0.2).

Set up the shooting method by hand, so you can control the convergence. The problem is that small changes in the initial conditions cause the solutions to blow up before r reaches R. This confuses FindRoot, which is being used to solve for the initial conditions.

There are two ways to proceed in V9; in V10 only the first method works. (1) Do it completely by hand with ParametricNDSolve or (2) use the "Shooting" method option "StartingInitialConditions". We get singularity warnings with both approaches when FindRoot happens to try the wrong initial conditions. One can use Quiet to suppress these messages, or live with them. I used Quiet below, but normally I just live with them.

### ParametricNDSolve (V10/V9)

sol = ParametricNDSolve[{eqs[[All, 1 ;; 2]], x1[R0] == a,
x2[R0] == b}, {x1, x2}, {r, R0, R}, {a, b},
"ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False}];

bcs = Quiet[
FindRoot[{x1[a, b][R] == Cb, x2[a, b][R] == Sp} /. sol, {a,
0.1}, {b, 12}],
ParametricNDSolve::ndsz
]
(*  {a -> 0.0000188874, b -> 20.8572}  *)

{x1i, x2i} = {x1[a, b], x2[a, b]} /. bcs /. sol;
{Cb, Sp}
{x1i[R], x2i[R]}
(*
{8, 30}
{8., 30.}
*)


### "StartingInitialConditions" (V9)

sol = Quiet[
NDSolve[eqs, {x1, x2}, {r, R0, R},
Method -> {"Shooting",
"StartingInitialConditions" -> {x1[R0] == 0.1, x1'[R0] == 0,
x2[R0] == 12, x2'[R0] == 0}}],
NDSolve::ndsz
];

{Cb, Sp}
{x1[R], x2[R]} /. sol
(*
{8, 30}
{{8.00087, 30.001}}
*)


### Update: "StartingInitialConditions" (V10)

The problem is that the solution {a -> 0.0000188807, b -> 20.8572} is close to the edge of the domain where solution is nonsingular for r < R = 0.4. The use of FindRoot in the built-in shooting method tends to step outside this domain if the starting value for x1[R0] == a is to large. In V10 (in this example at least), this seems to cause the shooting method to fail. The solution is to start with a value of a that is between the solution and the boundary of the domain. Instead of figuring all that out, it seems easier to use the first approach above, though.

sol = NDSolve[eqs, {x1, x2}, {r, R0, R},
Method -> {"Shooting",
"StartingInitialConditions" -> {x1[R0] == 0, x1'[R0] == 0,
x2[R0] == 20., x2'[R0] == 0}}]


The solution is (virtually) the same as the V9 solution.

• Did you use 9.0.1.0 for your second approach? It fails for me with 10.0.1.0 (with NDSolve::ndsz, which terminates the calculation) but not 9.0.1.0. Your first approach works flawlessly. Well done. – bbgodfrey Jan 29 '15 at 14:53
• @bbgodfrey Oops, you're right. I switched to V9.0.1 when I guessed that's what the OP was using. I will update to clarify. Thanks! – Michael E2 Jan 29 '15 at 17:20
• Do you plan to report this to Wolfram as a bug in 10.0.1.0? Certainly, that is how I would characterize it. – bbgodfrey Jan 29 '15 at 17:32
• Works great, thanks. – myguelangelo Jan 29 '15 at 18:28
• @bbgodfrey It certainly is a regression in this special case. But consider your own experience trying to solve this system, which seems similar to my own. (FindRoot above is rather sensitive to the starting values.) Given that BVPs in general are not guaranteed to have a solution, NDSolve must have some stopping criteria. Evidently they've changed. It may be that the change is an advance in 10 cases and a regression in this one. But you're right, it's probably worth reporting anyway. – Michael E2 Jan 29 '15 at 18:29

I think part of the issue lies with the coordinate singularity at r=0. I made the following transformations, first x1 -> x1/KO2 , x2->x2/KS . Then u1=x1/r and u2=x2/r. I also assumed that you were solving from R0>0 to R because you couldn't solve at R0=0. I set the boundary conditions u1[0]=0, u2[0]=0 . These transformations were used because

Laplacian[ f[r]/r, {r, θ, z}, "Spherical"] // Expand


(* = f''[r]/r *)

I found solutions for R=0.7

Then I wrote the equations as:

Clear[R, px, Mm, Ds, DO2, YXS, KS, Sp, Cb, KO2, R0, YXO, r, u1, u2, eqs2];
A = (Mm px)/(YXS DO2 KO2);
B = (Mm px)/(YXO Ds KS);
V = (r^3 u1[r])/(1 + r u1[r])  u2[r]/(1 + r  u2[r]);
Cbhat = Cb/KO2;
Sphat = Sp/KS;

eqs2 = {u1''[r] == A V, u1[0] == 0, u1[R] == R Cbhat, u2''[r] == B V,
u2[0] == 0, u2[R] == R Sphat}
R = .7; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.1; R0 = 0.000; YXO = 0.42857;

sol = NDSolve[eqs2, {u1, u2}, {r, 0, R}]
Plot[{(Evaluate[u1[r]] /. sol), R Cbhat}, {r, R0, R},  PlotRange -> {0, R Cbhat}]
Plot[{(Evaluate[u2[r]] /. sol), R Sphat}, {r, R0, R},  PlotRange -> {0, R Sphat}]


Beyond this I would try the shooting method suggested above. Beyond that I would try to scale u1, u2, and r to get rid of the very large parameter A. With the proper scaling a perturbative approach will suggest itself. If the solution actually ceases to exist for some value of R then an analytic approach will help. For starters, solutions disappear in pairs and at the bifurcation point the two solutions are equal. Thus, if you can find two solutions that look as if they're becoming equal as you vary R it would be indicative of a "saddle-node" bifurcation (of the sort that occurs at a=0 for the equation a-x^2=0).

Afterthoughts:

Clearly there is a boundary layer occurring for larger R. The fact that a numerical method fails to return a solution there is not surprising, nor is it necessarily indicative of a bifurcation point. There are various other tricks you might consider. This looks to be some sort of chemistry problem for which you're looking for stationary solutions? Either way you might consider solving the time-dependent version of the problem. It will always have a solution and if your steady solution is stable then the time dependent version will likely settle to it, assuming you have initial conditions that are in the basin of attraction for that solution.

COMMENT: MAJOR ERRORS CORRECTED IN THE BELOW and analytic solution broached

I should mention that you can reduce the dimension of the boundary value problem by noting that B u1"-A u2" = 0, so that w=B u1-A u2 obeys w" = 0 -> w = alpha r +beta and the boundary condtions give you beta=0 and alpha = B Cbhat - A Sphat. So you can keep the equation for u1 and replace u2= B/A u1 - (B Cbhat - A Sphat)/A . If you do this then you'll get the following equations:

Clear[R, px, Mm, Ds, DO2, YXS, KS, Sp, Cb, KO2, R0, YXO, r, u1, u2, eqs2];
A = (Mm px)/(YXS DO2 KO2)
B = (Mm px)/(YXO Ds KS);
V = (r^3 u1[r])/(1 + r u1[r]) u2[r]/(1 + r u2[r]);
V = V /. u2[r] -> ( B/A u1[r] - (B Cbhat - A Sphat) r/A)
(*
V = V/.u1[r]\[Rule]((Cbhat+Sphat)r+A/B u2[r])
*)
Cbhat = Cb/KO2;
Sphat = Sp/KS;

eqs2 = {u1''[r] == A V, u1[0] == 0, u1[R] == R Cbhat}
R = .8; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.1; R0 = 0.000; YXO = 0.42857;
A = 10^8
B
sol = NDSolve[eqs2, u1, {r, 0, R}];
u2[r_] = B/A ( Evaluate[u1[r]] /. sol) - (B Cbhat - A Sphat) r/A;
Plot[{(Evaluate[u1[r]] /. sol), R Cbhat}, {r, R0, R},  PlotRange -> {0, R Cbhat}, PlotStyle -> {Black, Red, Blue}]
Plot[{u2[r], Sphat r}, {r, 0, R}, PlotStyle -> {Black, Red}]


Note that I've overridden the defined value of A and set A=10^8. It wouldn't have mattered if I'd set it to 10^google. This is asymptopia for A. If you assume that u1 is bounded as A-> infinity (it is) then for u2 you get u2 = Sphat r in the large A limit. This still leaves you needing to solve:

u1" = A V but now, u2-> Sphat r in V. You'll have to figure out what the u2 solution means means but it means something in terms of the physics of your problem. If you pursue this you'd probably want to compute a correction for u2 as A=14,000 is not yet asymptotic. Now there is clearly a boundary layer and it is clearer but there's still work to be done to tease out the physics.

Below is the code that solves the 1-variable problem that remains after setting u2=Sphat r. I'll post this now but I think there are a few more easy steps to make.

Clear[R, px, Mm, Ds, DO2, YXS, KS, Sp, Cb, KO2, R0, YXO, r, u1, u2, eqs2];
A = (Mm px)/(YXS DO2 KO2)
B = (Mm px)/(YXO Ds KS);
V = (r^3 u1[r])/(1 + r u1[r]) u2[r]/(1 + r u2[r]);
V = V /. u2[r] -> ( Sphat r)

Cbhat = Cb/KO2;
Sphat = Sp/KS;

eqs2 = {u1''[r] == A V, u1[0] == 0, u1[R] == R Cbhat}
R = 1; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.1; R0 = 0.000; YXO = 0.42857;
A = 100000000
B
sol = NDSolve[eqs2, u1, {r, 0, R}];
u2[r_] = B/A ( Evaluate[u1[r]] /. sol) - (B Cbhat - A Sphat) r/A;
LogLogPlot[{(Evaluate[u1[r]] /. sol), R Cbhat}, {r, .00000001, R}, PlotRange -> {0, R Cbhat}, PlotStyle -> {Black, Red, Blue}]
Plot[{u2[r], Sphat r}, {r, 0, R}, PlotStyle -> {Black, Red}]

• I saw your solution and found it very interesting. Your assumption is roughly correct, it is a problem of diffusion within biofilms at steady state. I think the time-dependent version of the problem is somewhat more complicated, but that is what I'm working on now. I will take into consideration several of your advices, and I think I will need post more questions as I progress, thanks. – myguelangelo Jan 29 '15 at 19:34
• If you go to the time dependent problem then you'll find dw/dt = w"(r,t) where " = d^2/dr^2 . (all d's are partial derivatives). – JEP Jan 30 '15 at 1:22

Anther approach in V12.0 is to use the finite element method. For that we remove the zero derivatives (these are zero NeumannValues)

eqs = {{x1''[r] + 2/r x1'[r] == Vo/DO2,
x1[R] == Cb}, {x2''[r] + 2/r x2'[r] == Vs/Ds, x2[R] == Sp}}


with the remaining parameters:

Vs = px 1/YXS (Mm x2[r])/(KS + x2[r]) x1[r]/(KO2 + x1[r]);
Vo = px 1/YXO (Mm x2[r])/(KS + x2[r]) x1[r]/(KO2 + x1[r]);
R = 0.4; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.2; R0 = 0.000001; YXO = 0.42857;


We call NDSolveValue

sol = NDSolveValue[eqs, {x1, x2}, {r, R0, R},
Method -> "FiniteElement"];


Looking at the solution:

Through[sol[R]]
{8., 30.}
`

This what the other solutions found.