# solve tricky nonlinear ODE on a grid

I've been working on the same question for a while now and would appreciate any advice on what I'm doing wrong or any direction on how to actually solve the problem. :) On xzczd's advice, I retried with just using NDSolve. Also on xzczd's advice, I simplified the problem so that it is copy-paste-runnable.

Here's explicit code of what I'm doing.

(* Parameters *)
α = 4/100; γ = 2; ψ = 2; A = 1/10; θ = 10; δ = 0; σ = 1/10;

(* Define a function *)
ϕ[i_] := i - θ/2 i^2 - δ;

(* Boundaries *)
b = Rationalize[0.11619,10^-10]; (* Rationalize common-valued boundary point s*)

(* Some functions for ODE and solution to i *)
l1 = (1 - z) (1 - z y'[z]/y[z]);
l2 = z (1 + (1 - z) y'[z]/y[z]);
m = z^2 (1 - z)^2 y''[z]/y[z];

(* Solve for i as a function of y and y' *)
isolve = FullSimplify[Solve[{
(((A - i1) (1 - z) + (A - i2) z)/y[z])^(1/ψ) == α/(ϕ'[i1] (y[z] - z y'[z])),
(((A - i1) (1 - z) + (A - i2) z)/y[z])^(1/ψ) == α/(ϕ'[i2] (y[z] + (1 - z) y'[z]))}, {i1, i2}]];

(* There are three potential solutions here.  Only the first is
correct, I believe *)
Dimensions[isolve]
i1star = i1 /. isolve[[1, 1]];
i2star = i2 /. isolve[[1, 2]];
cstar = (A - i1star) (1 - z) + (A - i2star) z;

(* Simplify ODE *)
ode = Simplify[ ( α/(1 - 1/ψ) ((cstar/y[z])^(1 - 1/ψ) - 1) + ϕ[i1star] l1 + ϕ[i2star] l2 - γ /
2 (σ^2 l1^2 + σ^2 l2^2) + (σ^2 + σ^2)m/2) ] == 0;

(* Attempt at solution with NDSolve *)
sol = NDSolve[{ode, y[1/100] == b, y[99/100] == b}, y, {z, 1/100, 99/100}, Method -> {"Shooting", "StartingInitialConditions" -> {y[1/100] == b,
y'[1/100] == Rationalize[0.1536914, 10^-10]}}, WorkingPrecision -> 30]


As you can see, I'm starting to shoot from within the boundary, to sidestep the end-point singularities. The solution, I know, looks effectively (negative) quadratic, with endpoints at the boundaries. I was previously implementing a finite difference method because it seemed to deal with the end-points singularities better (I only solved on the interior of the grid). It also seemed to solve quickly.

My question is: is there a better, more robust way to solve this than shooting? Or is the math what it is and that's that?

Thank you.

• Well, you'd better add a specific working example rather than describe the problem in such an abstract way, or it's hard to give advice. With the information you've mentioned so far, I'll guess that the trouble might lie in FindRoot because solving a large nonlinear system with it can be troublesome and a good starting point is very important for FindRoot to give a proper answer. BTW, just to make sure, you need to solve the ODE yourself rather than using NDSolve, right? – xzczd Oct 28 '15 at 5:07
• Hi xzczd. Thanks for your reply. I was hoping someone like you was going to respond. I've read a lot of your online help about solving ODEs with NDSolve. I tried--A LOT-- with just NDSolve, but I found the shooting method to be unreliable in solving. I believe the problem is that the ODE is singular at both boundary points. By doing differencing myself, I was able to solve a simpler version of the problem and the differencing seemed to be much more robust. I'm still learning. – Jane Doe Oct 28 '15 at 19:05
• There's a typo in the definition of bi and I can't figure out how to fix it. And… currently your big code sample is a little frustrating, it'll be better if you can simplify the code sample a little, for example, the code calculating bi and b and ode isn't relevant, you can just give the result in the sample. If you're not sure whether your deductions for them are correct or not so want us to have a check (it's deprecated, of course), at least add a description for the original problem, I mean, you can express the problem with some equations (in traditional mathematical form). – xzczd Oct 29 '15 at 2:34
• This sounds like a related Q&A: mathematica.stackexchange.com/questions/91854/… – Michael E2 Oct 29 '15 at 14:23
• @xzczd It certainly looks complicated, but in part that's because it's not copy-pasteable. I don't have time to fix it up right now & explore. I think Jane Doe would be more likely to get help if the syntax errors were cleaned up. – Michael E2 Oct 30 '15 at 2:42

This computation is failing, because

(-((y[z] (y[z] + (1 - 2 z) y'[z])^2 (y[z] - z y'[z]))/(y[z] - (-1 + z) y'[z])^2))^(1/3)


is complex for typical initial conditions given in the question. Evidently, the wrong branch of cube root has been chosen. Instead, use Surd to assure that the real root is chosen (if it exists). To do so, modify the original ODE in the question as shown.

pode = ode /. Power[x__, 1/3] -> Surd[x, 3] /. Power[x__, 2/3] -> Surd[x, 3]^2
/. Power[x__, -1/3] -> 1/Surd[x, 3]
(* 1/100 (-8 + 4/5 2^(2/3) Sqrt[Surd[(y[z] (y[z] + (1 - 2 z) y'[z])^2 (y[z] -
z y'[z]))/(y[z] - (-1 + z) y'[z])^2, 3]/(y[z] (y[z] - z y'[z]))] - ...) *)


(The entire expression is quite long.) Then, after a bit of experimentation, the ODE with its boundary conditions can be solved by

sol = NDSolve[{pode, y[1/100] == b, y[99/100] == b}, {y, y'}, {z, 1/100, 99/100},
Method -> {"Shooting", "StartingInitialConditions" -> {y'[1/100] == .94777}}];
Plot[(y[z] /. sol), {z, 0.01, .99}, PlotRange -> All]


For completeness, the value y'[1/100] obtained by NDSolve is

(y'[z] /. sol) /. z -> .01
(* {0.947771} *)


It appears that a very good initial guess is needed, because the solution is near a separatrix of the ODE. The solution can be obtained in a more automated manner using approaches such as those given here.