7
$\begingroup$

I just got the new mathematica version 10 and tried to solve the following system of differential equation.

$$r^2\frac{d^2f}{dr^2} = 2f(1-f)(1-2f)+\frac{r^2}{4}h^2(f-1)$$ $$\frac{d}{dr}\left[r^2\frac{dh}{dr} \right]=2h(1-f)^2$$ $$ f(0)=h(0)=0, \quad f(\infty)=h(\infty)=1$$

and my code in mathematica,

Clear[eqn, bc, f, h]
eqn = {r^2*f''[r] == 2 f[r] (1 - f[r]) (1 - 2 f[r]) - r^2/4 (h[r])^2 (1 - f[r]),D[r^2*h'[r], r] == 
2 h[r] (1 - f[r])^2};
bc = {f[0.000001] == 0, f[20] == 1, h[0.000001] == 0, h[20] == 1};
sol = NDSolve[{eqn, bc}, {h, f}, {r, 0.000001, 20}]

Mathematica returned the error message:

NDSolve::ndsz: At r == 16.162462464292577`, step size is effectively zero; singularity or stiff system suspected. >>

Now, my problem is the New version 10 does not give me the InterpolatingFunction

while the previous version 8 gives me,

{{h->InterpolatingFunction[{{1.`*^-6,20.`}},"<>"],f->...}}

enter image description here

Does anybody know what's going wrong with the new version? How may I get the solution in V10? My sincere thanks for you patience and help!

It seems that the numerical solution behaves strangely in large $r$ region. If I set the boundary at $r=50$, the solution looks like this, enter image description here

Improve this question
$\endgroup$
8
  • $\begingroup$ can you show a plot of the solution you obtained for f(r) and h(r) in V8 for r=0..20? What options did you use with NDSolve for V8? Same as shown? or did you change something? $\endgroup$
    – Nasser
    Commented Jan 12, 2015 at 10:40
  • 1
    $\begingroup$ I get the warning in V9.0.1, but it still produces a solution: !Mathematica graphics $\endgroup$
    – Michael E2
    Commented Jan 12, 2015 at 11:17
  • $\begingroup$ Thank you @Michael E2 , that's exactly what I obtained from V8.0. But V10 gave me nothing $\endgroup$
    – mastrok
    Commented Jan 12, 2015 at 11:56
  • 1
    $\begingroup$ Maple solves this, using BVP method with mid-point, used it since there is singularity at t=0. Here is screen shot. !Mathematica graphics Tried all the mid points methods in M with NDSolve, but none of them worked. same error. V 10.02 $\endgroup$
    – Nasser
    Commented Jan 12, 2015 at 12:26
  • $\begingroup$ The output I got from V8 is the interpolatingFunction which I can use them {{h->InterpolatingFunction[{{1.*^-6,20.}},"<>"],f->...}} when I plot them out I just got the same plot as @Michael E2. However V10 does not give me that $\endgroup$
    – mastrok
    Commented Jan 12, 2015 at 12:29

1 Answer 1

Reset to default
5
$\begingroup$

The problem is that NDSolve is using the shooting method and for some of the initial conditions it tests lead to a singularity. Indeed there seems to be a small neighborhood around the desired solution where the integration does not blow up before r == 20. I don't know why previous versions of Mathematica were able to push through and find an adequate solution. In any case V10 seems to give up when it gets an error.

One approach is to transform the differential equation to be a bit better behaved near r == 0. We can replace f[r] by r^2 g[r] and h[r] by r j[r].

Clear[eqn, bc, f, h, g, j]
With[{r1 = 0.000001, r2 = 20},
 eqn = {r^2*f''[r] == 2 f[r] (1 - f[r]) (1 - 2 f[r]) - r^2/4 (h[r])^2 (1 - f[r]), 
   D[r^2*h'[r], r] == 2 h[r] (1 - f[r])^2};
 bc = {f[r1] == 0, f[r2] == 1, h[r1] == 0, h[r2] == 1};
 {sol} = NDSolve[{eqn, bc} /. {f -> (#^2 g[#] &), 
     h -> (# j[#] &)}, {g, j}, {r, r1, r2}]
 ]

Plot[{r j[r], r^2 g[r]} /. sol // Evaluate, {r, 1.*^-6, 20}]

Mathematica graphics

Improve this answer
$\endgroup$
1
  • $\begingroup$ In fact I have a second, it seems that there is some problem with the numerical solution if I set $r_2$ to a large value, say, $r_2 =50$. The numerical solution behaves strangely at large $r$ region and the shape at small $r$ also change. I will upload the screenshot in my question. $\endgroup$
    – mastrok
    Commented Jan 13, 2015 at 5:33

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.