# Why does NDSolve blow up when given my ODE but bvp4c in Matlab does not?

I am numerically solving the following ODE initially using NDSolve in Mathematica(updated and corrected):

$$-(-\frac{z'(r)}{r \sqrt{z'(r)^2+1}}-\frac{z''(r)}{\left(z'(r)^2+1\right)^{3/2}})=A_1(z(r)+H)+A_2(\frac{A_3}{\sqrt{z(r)^2+r^2}-1})^3$$

subject to boundary conditions of $$z'(0) = 0$$, and $$z(\infty) = - H$$. In this case, $$r = 4$$ or even 1 can be treated as infinity.

The left-hand-side is twice the mean curvature of a surface of revolution characterized by $$z(r)$$ in cylindrical coordinates. Parameters $$A_1, A_2, A_3$$ and $$H$$ are provided. I anticipated a solution that looked like this:

The expected solution (blue curve) is perturbed by the presence of the sphere near $$r = 0$$ and resumes its unperturbed shape (a straight line) at $$r = \infty$$. Note that H is the vertical distance of z(r) to the line of $$z = 0$$ at $$r = \infty$$..

The Mathematica code used:

A1 = 200;
A2 = 1.86*10^7;
A3 = 0.0002;
H=1;
f[z_, r_] := A1 (z + H) + A2 (A3/(Sqrt[z^2 + r^2]-1))^3;
k := -(z''[r]/(z'[r]^2 + 1)^(3/2)) - If[r == 0, 0, z'[r]/(r Sqrt[z'[r]^2 +1])]);
sol = NDSolve[{-k == f[z[r], r], z'[0] == 0, z[1] == -H}, z, {r, 0, 1}];


However, NDSolve was not able to yield sensible results, with $$z(r)$$ either blowing up to $$z = \infty$$ or simply "encountered stiffness", even if variations of the boundary conditions were attempted, i.e.

1) z'[0.0001] = 0, z[0.0001] = -1.01

2) z'[1] = 0, z[1] = - H

3) z'[0.0001 = 0, z[1] = -H

To my surprise, Matlab on the other hand handled the identical system well using bvp4c (4th order RK method) without blow-up, and yielded the solution shown in the figure above.

Any clue as to why Mathematica resulted in a blow-up solution yet Matlab converged well? Any explanations will be greatly appreciated.

• Could you show your Matlab code also? – GerardF123 Nov 8 '18 at 15:25
• Sorry I'm new to StackExchange, how do I send you my Matlab code? – Linmin Nov 8 '18 at 19:37
• You should be able to paste it here I would think. – GerardF123 Nov 8 '18 at 20:13
• clc xlow=1e-10; xhigh=2; solinit = bvpinit(linspace(xlow,xhigh,20),[-1.03 0]); sol = bvp4c(@twoode,@twobc,solinit); xint = linspace(xlow,xhigh); Sxint = deval(sol,xint); plot(xint,Sxint(1,:),xint,-(1-xint.^2).^0.5) pbaspect([1 1 1]) xlabel('r') ylabel('z') %% --- function dydx = twoode(x,y) A0=5.36e-14; Lambda=2.0e-10; G=1.0e4; R=1.0e-6; Tm=933; H=1e-6; A1=GR^2/A0/Tm; A2=GH*R/A0/Tm; A3=R/A0; A4=Lambda/R; % A1=0; % A2=0; % A3=0; % A4=0; zpp=(A1*y(1)+A2+A3*(A4/(sqrt(y(1)^2+x^2)-1))^3)*(1+y(2)^2)^1.5-y(2)*(1+y(2)^2)/x; %zpp=-abs(y(1)); dydx = [ y(2) zpp ]; end – Linmin Nov 8 '18 at 21:17
• %% ---------------This is the second half------------------------ function res = twobc(ya,yb) H=1e-6; R=1.0e-6; res = [ ya(2) yb(1)+H/R]; end – Linmin Nov 8 '18 at 21:18

First, there is a mistake in the sign k. And secondly, for such tasks we use a special method that expands the possibilities of the shooting method. I will demonstrate the working code

A1 = 200;
A2 = 186*10^5;
A3 = 1/5000;
a = Rationalize[A2*A3^3];
H = 1; r0 = 10^-5;
f[z_, r_] := A1*(z + H) + a*(1/(Sqrt[z^2 + r^2] - 1))^3;
k = -z''[r]/(z'[r]^2 + 1)^(3/2) - z'[r]/(r Sqrt[z'[r]^2 + 1]);
sol = ParametricNDSolveValue[{k == f[z[r], r], z'[r0] == 0,
z[r0] == z0}, z, {r, r0, 1}, {z0}];
Plot[Evaluate[Table[sol[z0][r], {z0, -1.15, -1.05, .01}]], {r, r0, 1},
PlotRange -> All]


Find the parameter that satisfies the boundary condition on the right border

FindRoot[sol[z0][1] == -H, {z0, -1.05}]

Out[]= {z0 -> -1.03176}
{Plot[sol[-1.03176][r], {r, r0, 1}],Plot[{sol[-1.03176][r], -Sqrt[1 - r^2]}, {r, r0, 1}]}


The author insists that the correct model corresponds to his code. We give a solution for this case. The solution method does not change.

A1 = 200;
A2 = 186*10^5;
A3 = 1/5000;
a = Rationalize[A2*A3^3];
H = 1; r0 = 10^-5;
f[z_, r_] := A1*(z + H) + a*(1/(Sqrt[z^2 + r^2] - 1))^3;
k = z''[r]/(z'[r]^2 + 1)^(3/2) + z'[r]/(r Sqrt[z'[r]^2 + 1]);
sol = ParametricNDSolveValue[{k == f[z[r], r], z'[1] == z1,
z[1] == -H}, z, {r, r0, 1}, {z1}, WorkingPrecision -> 30];
{Plot[Evaluate[
Table[sol[z1][
r], {z1, .00024816068, .00024816069, .000000000001}]], {r, r0,
1}, PlotRange -> All],
Plot[Evaluate[
Table[sol[z1][
r], {z1, .00024816068, .00024816069, .000000000001}]], {r,
r0, .01}, PlotRange -> All]} // Quiet


• Thanks Alex. I have double-checked the literature that the sign of k in my original post is correct, it is indeed -k = f (z(r) ,r ), and the signs for z'' and z' in the expression of k are both negative, but the wavy behavior in your post is not expected in the actual solution. What advantage does ParametricNDSolveValue have over just NDSolve? I see that NDSolve could equally handle well k = f (z,r), and neither of them can solve - k = f(z,r) without blowing up. – Linmin Nov 7 '18 at 15:11
• Maybe you should check three times. Look at the original equation i.stack.imgur.com/yP3xE.gif and your code. What solution do you expect to get? – Alex Trounev Nov 7 '18 at 15:50
• Opps you are correct, the expression in the image is incorrect. The expression in the code is what I am supposed to solve, and is what lead to blow-up issue. – Linmin Nov 7 '18 at 16:11
• Are you kidding? Give a link to the article or book where you took this equation. – Alex Trounev Nov 7 '18 at 16:18
• The equation I took is from this link. link Eqn 47 and 48, the definitions of parameters are given directly below Eqn 48, but you may use the values I provided initially. The solution z(r) is shown in Fig 3 (b) via 4th order RK method. I can provide a PDF copy if you don't have access to Journal of Computational Physics. Thanks! – Linmin Nov 7 '18 at 16:37

This type of problem has been discussed in this site for many times, for example this, and you can find more by searching Shooting in this site. For completness I'd like to add a solution based on finite difference method (FDM). I'll use pdetoae for the task:

A1 = 200;
A2 = 186/100 10^7;
A3 = 2 10^-4;
H = 1;
f[z_, r_] = A1 (z + H) + A2 (A3/(Sqrt[z^2 + r^2] - 1))^3;
k = -(z''[r]/(z'[r]^2 + 1)^(3/2)) - z'[r]/(r Sqrt[z'[r]^2 + 1]);

{eq, bc} = {-k == f[z[r], r], {z'[0] == 0, z[1] == -H}};

(* Remove the removable singularity *)
neweq = 0 == (Together[Subtract @@ eq] // Numerator);

points = 25; domain = {0, 1}; grid = Array[# &, points, domain]; difforder = 4;
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[z[r], grid, difforder];
del = #[[2 ;; -2]] &;
ae = ptoafunc@neweq // del;
aebc = ptoafunc@bc;
initial[r_] = -1;
sollst = FindRoot[{ae, aebc}, Table[{z[r], initial@r}, {r, grid}],
WorkingPrecision -> 16][[All, -1]]
sol = ListInterpolation[sollst, grid]

Plot[sol[r], {r, 0, H}]


• Thanks very much, I will be reading those posts about shooting method for this peculiar ODE system. I also noticed that you and the first answer rationalized all parameters even if I originally assigned 0.0002 to A3, which seemed counterintuitive to me since my common sense told me that float-point numbers would speed up numerical calculations. Any reason why? – Linmin Nov 8 '18 at 21:26
• @Linmin High level numeric functions like FindRoot, NDSolve automatically turns to float-point numbers by default, so you don't need to manually change the parameters to float-point number. On the other hand, using float-point number obstructs adjustion for options like WorkingPrecision. (Just remove this option from my code and use A3=0.0002 and retry. ) BTW, you can also use Rationalize[…, 0] for rationalizing. – xzczd Nov 9 '18 at 3:23
• Thanks a lot, that does help clarify some points. – Linmin Nov 9 '18 at 22:55