# Solving non-linear differential equation

For my research project I need to solve the following integral

I =
1/2 Integrate[π Sqrt[r[θ]^2 Sin[θ]^2]
Sqrt[r[θ]^2 + Derivative[1][r][θ]^2/
(1 + 1/r[θ]^2 - 3/r[θ] + r[θ]^2)],
{θ, 0, 0.005}]


where r[θ] is a solution of following differential equation

-2 r[θ] + 12 r[θ]^2 - 22 r[θ]^3 + 12 r[θ]^4 - 6 r[θ]^5 + 12 r[θ]^6 -
4 r[θ]^7 - 2 r[θ]^9 +
(Cot[θ] r[θ]^2 - 3 Cot[θ] r[θ]^3 + Cot[θ] r[θ]^4 + Cot[θ] r[θ]^6)
Derivative[1][r][θ] +
(-2 r[θ] + (15 r[θ]^2)/2 - 3 r[θ]^3 - 4 r[θ]^5) Derivative[1][r][θ]^2 +
Cot[θ] r[θ]^2 Derivative[1][r][θ]^3 +
(r[θ]^2 - 3 r[θ]^3 + r[θ]^4 + r[θ]^6) Derivative[2][r][θ]


with initial conditions {r[0] == r0, r'[0] == 0}.

It's a second order non-linear differential equation. The independent variable θ runs from θ = {0, 0.005}. There is a singularity at θ = 0, which can be avoided by taking θ = 10^-10 (some small value). Then initial conditions modify to {r[10^-10] == r0, r'[10^-10] == 0}.

I have to choose r0 such that by solving the above equation one should get r[0.005] = 10000. The initial guess for r0 (depending upon the working precision) can be around r0 = 199.958344.

At the end, integral should gives value I = 76.96884.

I am solving this problem in Mathematica and getting answer like I = 76.9844.

• Your differential equation isn't actual an equation. Equations should have == somewhere. – Sjoerd C. de Vries Oct 9 '15 at 10:25
• Where you write (r^′′)[θ] in the above, do you really mean Derivative[2][r][θ]? – m_goldberg Oct 9 '15 at 10:40
• I'm set WorkingPrecision -> 30 and Integral is: I = 76.9843992958913645452412310988. Maybe this is the correct answer. – Mariusz Iwaniuk Oct 9 '15 at 10:52
• Dear Sjoerd C. de Vries, I have just given the expression for the Left hand side. The right hand side of the equation is zero. Thanks for pointing out. – chandra Oct 9 '15 at 11:14
• Dear @m_goldberg, Thank you for your time. Can you please post the full mathematica code here. I have also used WorkingPrecision -> 30, but didn't get the answer. Thanks in advance. – chandra Oct 9 '15 at 11:18

The following leads me to believe that the value of the integral as given by I_Mariusz is correct.

First, solve the ODE itself with MaxStepFraction -> 10^-6.

ϵ = 10^-50;
{s, t} = NDSolveValue[{-2 r[θ] + 12 r[θ]^2 - 22 r[θ]^3 + 12 r[θ]^4 - 6 r[θ]^5 +
12 r[θ]^6 - 4 r[θ]^7 - 2 r[θ]^9 + (Cot[θ] r[θ]^2 - 3 Cot[θ] r[θ]^3 + Cot[θ] r[θ]^4 +
Cot[θ] r[θ]^6) Derivative[1][r][θ] + (-2 r[θ] + (15 r[θ]^2)/2 - 3 r[θ]^3 -
4 r[θ]^5) Derivative[1][r][θ]^2 + Cot[θ] r[θ]^2 Derivative[1][r][θ]^3 + (r[θ]^2 -
3 r[θ]^3 + r[θ]^4 + r[θ]^6) Derivative[2][r][θ] == 0, r'[ϵ] == 0, r[5/1000] == 10000},
{r, r'}, {θ, ϵ, 5/1000}, WorkingPrecision -> 50, MaxSteps -> 10^7,
MaxStepFraction -> 10^-6, Method -> {"Shooting",
"StartingInitialConditions" -> {r'[ϵ] == 0, r[ϵ] == 199958364/1000000}}];


which yields (after a few hours of calculation) the interpolation functions for r and r'.

LogPlot[{s[θ], t[θ]}, {θ, ϵ, .005}, PlotRange -> All]


Note that the solution is growing faster than exponentially near θ -> 0.005. The initial (called r0 in the question) and final values of r are

{s[ϵ], s[5/1000]}
(* {199.95836405695030604588498996975119999174071797848,
10000.000000000000000005072133417578396839834109} *)


To estimate the accuracy of this solution, repeat the computation with MaxStepFraction -> 10^-5, which takes only several minutes, and compare the resulting solution, designated {s5, t5}, to {s, t}.

LogPlot[{Abs[s[θ] - s5[θ]]/s5[θ], Abs[t[θ] - t5[θ]]/t5[θ]}, {θ, ϵ, .005}, PlotRange -> All]


Agreement is excellent. The integral can be computed by

π/2 NIntegrate[s[θ] Sin[θ] Sqrt[s[θ]^2 + t[θ]^2/(1 + 1/s[θ]^2 - 3/s[θ] + s[θ]^2)],
{θ, ϵ, 0.005}]
(* 76.98439929825058 *)


or by

sx = Flatten[s["Grid"]]; sy = s["ValuesOnGrid"]; ty = t["ValuesOnGrid"];
int = Interpolation[Table[{sx[[i]], sy[[i]] Sin[sx[[i]]] Sqrt[sy[[i]]^2 +
ty[[i]]^2/(1 + 1/sy[[i]]^2 - 3/sy[[i]] + sy[[i]]^2)]}, {i, Length[sx]}]];
π/2 Integrate[int[θ], {θ, ϵ, 0.005}]
(* 76.98439929486328 *)


The two results agree to ten significant figures. I also have used NIntegrate with WorkingPrecision -> 50, MaxRecursion -> 50, as well as with other values of these two quantities. Although NIntegrate then generates error messages, it still gives results agreeing with those above to nine significant figures.

I believe that we may safely conclude that the answer is approximately 76.9843992.

• "after a few hours of calculation" - wow, that hard, eh? – J. M. is away Oct 10 '15 at 15:39
• @bbgodfrey, Thank you very much for your help. I think what you are saying is correct. And now I am convince that the numerical value of the integral should be I=76.98439 instead of I = 76.96884. Thank you – chandra Oct 10 '15 at 15:59

At first we find r0:

 eq = {-2 r[\[Theta]] + 12 r[\[Theta]]^2 - 22 r[\[Theta]]^3 +
12 r[\[Theta]]^4 - 6 r[\[Theta]]^5 + 12 r[\[Theta]]^6 -
4 r[\[Theta]]^7 - 2 r[\[Theta]]^9 + (Cot[\[Theta]] r[\[Theta]]^2 -
3 Cot[\[Theta]] r[\[Theta]]^3 + Cot[\[Theta]] r[\[Theta]]^4 +
Cot[\[Theta]] r[\[Theta]]^6)*r'[\[Theta]] + (-2 r[\[Theta]] + (15 r[\[Theta]]^2)/2 -
3 r[\[Theta]]^3 - 4 r[\[Theta]]^5)*(r'[\[Theta]])^2 +
Cot[\[Theta]] r[\[Theta]]^2*(r'[\[Theta]])^3 + (r[\[Theta]]^2 -
3 r[\[Theta]]^3 + r[\[Theta]]^4 + r[\[Theta]]^6)*(r''[\[Theta]]) == 0}

sol2 = With[{\[Epsilon] = 10^-50},
r /. ParametricNDSolve[{eq, r[\[Epsilon]] == r0,
r'[\[Epsilon]] == 0}, r, {\[Theta], \[Epsilon], 5/1000}, {r0},
WorkingPrecision -> 50, MaxStepFraction -> 0.00001,
MaxSteps -> 10^6]];
sol3[r0_?NumericQ] := sol2[r0][5/1000];
SetPrecision[
FindRoot[sol3[r0] == 10000, {r0, 199.957, 199.959}, Method -> "Brent",
WorkingPrecision -> 50], 50]


$\{r0\to 199.95836405695030604588494221489902476795675126277\}$

With Maple I have got r0 = 199.96131889676149,a very big difference in results.

Calculate the integral:

 sol = First@With[{\[Epsilon] = 10^-30},r /. NDSolve[{eq,
r[\[Epsilon]] == 199.95836405695030604588494221489902476795675126276935236196857881508093072408836, r'[\[Epsilon]] == 0},
r, {\[Theta], \[Epsilon], 5/1000}, MaxStepFraction -> 0.00001,
MaxSteps -> 10^6, WorkingPrecision -> 30]];
\[Epsilon] = 10^-30;
i = (1/2)*NIntegrate[Pi*Sqrt[sol[\[Theta]]^2 Sin[\[Theta]]^2]*Sqrt[sol[\[Theta]]^2 +
sol'[\[Theta]]^2/(1 + 1/sol[\[Theta]]^2 - 3/sol[\[Theta]] +
sol[\[Theta]]^2)], {\[Theta], \[Epsilon], 5/1000},WorkingPrecision -> 30, MaxRecursion -> 10000]


$i = 76.9843992958929803977886284771$

Plot:

Plot[{sol[\[Theta]], 199.95836}, {\[Theta], 0, 0.0051},
Prolog -> {Line[{{0.005, 0}, {0.005, 2000}}]},
PlotRange -> {{0, 0.0051}, {0, 2000}},
PlotLegends -> {"r[\[Theta]]", "r0"}]


With Maple I have got i = 76.9843992970933.In both programs results are almost similar.

• What do you mean by your last sentence? Also, how did you obtain I = 76.96884 in your question and how do you know that it is the right answer? – bbgodfrey Oct 9 '15 at 19:49
• @Dear I_Mariusz, Thank you very much for trying. I have wrote the same kind of code which you have posted here and getting the same answer for Integral. It seems to me that the problem is with calculating r0. Unfortunately the published version of the paper, which I am trying to reproduce, has I = 76.96884. However, authors in that paper using matlab's bvp4c. Since, we are using Mathematica, may be the problem lies here. If you have any comments in this regard, I will highly appreciate that. – chandra Oct 10 '15 at 3:10
• Dear bbgodfrey, Please check the mathematica file posted by I_Mariusz to get r0 and I. For the answer, I am following a published paper where, I = 76.96884 was archived. – chandra Oct 10 '15 at 3:14
• @Chandra Can you provide a link to the paper as it all sounds confusing as to whether I=76.96884 or not. What is correct and what is wrong. U said it is "unfortunate"...why? Sometimes ambiguous input can waste precious time on part of those who are trying to help. – thils Oct 10 '15 at 3:59
• @bbgodfrey.I do not know whether this is the correct answer.My last sentence->We can skip and it is not no big deal.Sorry my English bad.(Google Translator not perfect yet) – Mariusz Iwaniuk Oct 10 '15 at 6:40