3
$\begingroup$

I'm trying solve this problem:

g'(r) = a(r)g(r)/r, (1/r)a'(r) = g(r)^2-1

which have the following boundary conditions: a(0)=n and g(infinity)=1, where n>=1.

I have problem with the error for n=1[hence, failed for the other n values as well.]:

"Infinite expression 1/0.^2 encountered"

using DSolve directly.

A work colleague solved for the values n=1,2,5,10 using the software MAPLE. But, I don't know to use MAPLE and I would like to get these results by Mathematica. Below are the plots found by MAPLE.enter image description here

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ Please include the Mathematica code text in your question. $\endgroup$
    – xzczd ♦
    Commented Aug 8, 2018 at 16:58
  • 3
    $\begingroup$ Welcome to Mma.SE. Start by taking the tour now and learning about asking and what's on-topic. Always edit if improvable, show due diligence, give brief context, include minimal working example of code and data you are working on in formatted form. By doing all this you help us to help you and likely you will inspire great answers. The site depends on participation, as you receive give back: vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned $\endgroup$
    – rhermans
    Commented Aug 8, 2018 at 17:10
  • 1
    $\begingroup$ See this previous answer. Compactifying in MMA is rather trivial. $\endgroup$
    – Hector
    Commented Aug 8, 2018 at 17:10
  • 1
    $\begingroup$ I find that DSolve does not produce an error message but instead returns unevaluated, meaning that it cannot solve the problem.. $\endgroup$
    – bbgodfrey
    Commented Aug 8, 2018 at 21:34
  • 1
    $\begingroup$ Maple code:sol := proc (n) options operator, arrow; dsolve([(D(g))(r) = a(r)*g(r)/r, (D(a))(r)/r = g(r)^2-1, a(1/100000) = n, g(10) = 1], numeric) end proc;plots:-display(seq({plots:-odeplot(sol(n), [r, g(r)], linestyle = n, color = ColorTools:-Color([1/n, 1/(1+n), 1/(n+4)]))}, n = 1 .. 7));plots:-display(seq({plots:-odeplot(sol(n), [r, a(r)], linestyle = n, color = ColorTools:-Color([1/n, 1/(1+n), 1/(n+4)]))}, n = 1 .. 7)); $\endgroup$
    – Mariusz Iwaniuk
    Commented Aug 11, 2018 at 17:17

2 Answers 2

Reset to default
4
$\begingroup$

Another way:

eps = 10^-10; end = 10;
sol[n_] := NDSolve[
{g'[r] == a[r]*g[r]/r, a'[r]/r == g[r]^2 - 1, a[eps] == n, g[end] == 1}
,{a[r], g[r]}, {r, eps, end}, Method -> {"Shooting", 
"StartingInitialConditions" -> {a[end] == eps, g[end] == (1 - eps)}}];

Plot[Evaluate[g[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end},PlotRange -> All, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "g(r)"]
Plot[Evaluate[a[r] /. sol[#] & /@ {1, 2, 5, 10}], {r, eps, end}, PlotRange -> All, 
PlotLegends -> {"n=1", "n=2", "n=5", "n=10"}, PlotLabels -> "a(r)"]

enter image description here

For larger values e.g.end=20:

eps = 10^-20; end = 20;
sol[n_] := NDSolve["Code the same",WorkingPrecision -> 25];
Plot["Code the same"] // Quiet(*Error messages form InterpolatingFunction to Quiet *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for your answer, you've been a big help. $\endgroup$
    – Andrey
    Commented Aug 12, 2018 at 18:49
  • $\begingroup$ Dear @MariuszIwaniuk, I was very grateful for the help on this occasion. Now, if it is possible to help again, how could I make a projection of these solutions in the plan? That is, would it be possible to make a density plot? $\endgroup$
    – Andrey
    Commented Jul 19, 2019 at 19:19
4
$\begingroup$

Error in boundary condition corrected

To solve this ODE system, begin by determining its behavior near r == 0.

Series[{g'[r] - (a[r] g[r])/r, a'[r]/r + 1 - g[r]^2}, {r, 0, 2}] // Normal;
Thread[Flatten@CoefficientList[r % /. a[0] -> n, r] == 0];
Solve[%, {a'[0], a''[0], a'''[0], a''''[0], g[0], g'[0], g''[0], g'''[0]}]

(* a'[0] -> 0, a''[0] -> -1, a'''[0] -> 0, a''''[0] -> 0,
   g[0] -> 0, g'[0] -> 0, g''[0] -> 0, g'''[0] -> 0  *)

In other words, g[r] vanishes to all orders at the origin for most values of n. However, if n is a positive integer, the same set of commands leaves Derivative[n][g][0] undetermined. So, g[eps] is approximately proportional to eps^n. Next, solve the equations numerically by shooting, using the method previously employed here and elsewhere on this site. (As noted in my comment above, DSolve returns unevaluated when applied to these equations.) A typical calculation for n == 5 is given by

eps = 10^-5; end = 20;
sp = ParametricNDSolveValue[{g'[r] == a[r] g[r]/r, a'[r]/r == g[r]^2 - 1, 
    a[eps] == n0, g[eps] == eps^n0 gp, 
    WhenEvent[g[r] > 101/100, {bool = 1, "StopIntegration"}], 
    WhenEvent[a[r] < 0, {bool = 0, "StopIntegration"}]}, 
    {a[r], g[r]}, {r, eps, end + 5}, {gp, wp0, n0}, 
    WorkingPrecision -> wp0, Method -> "StiffnessSwitching"];
bl = 0; bu = 2; imax = 100; wp = 75; n = 5;
Do[bool = -1; bmiddle = (bl + bu)/2; st = sp[bmiddle, wp, n]; 
    rm = (First[st] /. r -> "Domain")[[1, 2]]; 
    If[bool == 0, bl = bmiddle, bu = bmiddle]; ip = i; 
    If[bool == -1, Return[]], {i, imax}]; 
s[n] = st;

For larger values of n, wp must be increased. For instance, with n = 10, wp = 90 is required. Results for n = {1, 2, 5, 10} are plotted below.

Plot[{s[1], s[2], s[5], s[10]}, {r, eps, 10}, PlotRange -> All, AxesLabel -> {r, "a, g"}, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, Medium}]

enter image description here

The solutions takes no more than several seconds each.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I have appreciated your answer. Thank you very much @bbgodfrey. You've already been very helpful. A friend got solve also for values n>1, but he used the software MAPLE. I don't know to use the MAPLE, but I'll try to use your answer like a procedure. $\endgroup$
    – Andrey
    Commented Aug 10, 2018 at 1:49
  • $\begingroup$ @AndréCavalcante Because the ODEs are nonlinear, not all boundary conditions are consistent with them. In this case, I believe that A[0] == 1 is the only consistent boundary condition, unless you are willing to accept g[r] == 0 as a solution. Perhaps, I could get Mathematica to produce numerical solutions for other values of n, but they would be incorrect, I believe. $\endgroup$
    – bbgodfrey
    Commented Aug 10, 2018 at 1:58
  • $\begingroup$ I have edited my question and added the plots obtained by MAPLE. Indeed, using its answer it was not possible to solve for values greater than n == 1. In addition, by behavior near r == 0, I have noticed that for n == 2, for example, leaves g''[0] undetermined. $\endgroup$
    – Andrey
    Commented Aug 11, 2018 at 16:58
  • $\begingroup$ @AndréCavalcante I just realized that I had made an error in the g[eps] boundary condition. In fact, any positive integer value of n gives a valid solution. Sorry. $\endgroup$
    – bbgodfrey
    Commented Aug 11, 2018 at 18:17
  • $\begingroup$ Again, thanks for your help. $\endgroup$
    – Andrey
    Commented Aug 12, 2018 at 18:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.