Finding out the nodes from spherical symmetric equation for mathematica and then plotting it in gnu-plot

439 Views Asked by At

What I want to do is to draw a curve for the spherical partial differential equation for S at rho=0:

\begin{equation} \frac{\partial^2S}{\partial \rho^2}+\frac{d-1}{\rho}\,\frac{\partial S}{\partial \rho} -S+S^3=0 \end{equation}

I have drawn by mathematica by applying the code:

    Manipulate[
 Plot[Evaluate[
   S[r] /. (NDSolve[{S''[r] + 1/r S'[r] - S[r] + S[r]^3 == 0, 
         S[1/1000] == #, S'[1/1000] == 0}, S, {r, 1/100, 15}] & /@ 
      Range[ic - wi, ic + wi, wi/2])], {r, 1/100, 10}, 
  PlotRange -> {-2.5, 5}, Frame -> True, AspectRatio -> 1/2, 
  ImageSize -> 800], {{ic, 3.5, "bundle start"}, -7, 
  7}, {{wi, 1.5, "bundle width"}, .5, 2}]

the code give the graph like: enter image description here

But the graph will be according to the article: enter image description here

The professor mailed me to do in this way but I'm bit confuse what exactly he says.

The professor suggests me this

The nodes of a curve are where its values cross zero. So the values of rho where S=0 in this case. At the mathematica file you should plot only one curve. When the central value, S at rho=0, is slightly below 2.2062 for 2+1 dimensions, then going outwards the curve approaches zero, but then diverges to plus infinity. If you are slightly below that value, it diverges to minus infinity. Fine tuning, i.e. shooting, between this two behavior there is one solution where the value tends to zero at infinity. If you do the same around the central value 3.3319, then you get a solution which crosses the S=0 line before it tends to zero, and then blows up at the positive or negative direction. The limiting solution has one node at around rho=3.7 and tends to zero at infinity.

After he said to use gnuplot.

So how do I get the the Fig with different nodes and the value of TABLE:1 (in the paper) with the Figure (3) and (4)

I have tried to get node for single curve by using the mathematica code but the central values given in the TABLE:1 page no 8 (below the page) aren't resembling with mine. I'm I doing wrong to get the nodes?

enter image description here

Source: article page no:7, Fig: no 1,2

1

There are 1 best solutions below

3
On

The boundary conditions are important. Your code has (approximate) conditions on $S'(0)$ and you vary the value at $S(0)$ You have to find some way to enforce that $S\rightarrow 0$ as $\rho\rightarrow\infty$. Your professor is giving you one way to do that. He is suggesting that you search for the right values of $S(0)$ and $S'(0)$ that give you solutions that decay at infinity. The authors of the paper searched through the right values and tabulated some of them for $S(0)$ to 8 digits of precision. They don't give you their results, however, for $S'(0)$.

It is an optimization problem, but it's basically "trial and error" done in a systematic way. The authors of the paper have done that for you for this equation. They found the stable (that don't blow up) solutions to the differential equation, and tell you that there are a family of stable solutions, each of higher and higher oscillatory order. They call this oscillation number as the number of "nodes".

If you insist on verifying the results of that paper, you'll have to write your own optimization routine that find the right $S'(0)$ and $S(0)$ values that cause $S\rightarrow 0$ as $\rho\rightarrow 0$. You'll get the numbers in the paper if you do this. I'm sorry I can't be more specific, because this is a non-trivial coding task. But it can be done.

Edit

The below gives you one plot from the paper. I won't be editing this post any further about gnuplot, or how to extend this to $D=3$, or how to FIND the initial conditions that satisfy the "vanish at infinity" boundary condition. I described how you do that above, it is simply a systematic trial and error, just as your professor described. Best of luck.

S0 = S[r] /. 
   NDSolve[{S''[r] + 1/r S'[r] - S[r] + S[r]^3 == 0, 
     S[1/10000000000000000] == 2.20620086, 
     S'[1/10000000000000000] == 1}, S, {r, 1/1000, 15}];
S1 = S[r] /. 
   NDSolve[{S''[r] + 1/r S'[r] - S[r] + S[r]^3 == 0, 
     S[1/10000000000000000] == 3.33198927, 
     S'[1/10000000000000000] == 0}, S, {r, 1/1000, 15}];
S2 = S[r] /. 
   NDSolve[{S''[r] + 1/r S'[r] - S[r] + S[r]^3 == 0, 
     S[1/10000000000000000] == 4.15009404, 
     S'[1/10000000000000000] == 0}, S, {r, 1/1000, 15}];
S3 = S[r] /. 
   NDSolve[{S''[r] + 1/r S'[r] - S[r] + S[r]^3 == 0, 
     S[1/10000000000000000] == 4.82960282, 
     S'[1/10000000000000000] == 0}, S, {r, 1/1000, 15}];
Plot[{S0, S1, S2, S3}, {r, 1/100, 14}, PlotRange -> Full]