3
$\begingroup$

Okay, I have this not so pretty 2nd order non-linear ODE I should be able to solve numerically.

$$f''(R) + \frac{2}{R} f'(R)=\frac{0.7}{R} \left( \frac{1}{\sqrt{f(R)}} - \frac{0.3}{\sqrt{1-f(R)}} \right),$$

$$f(1)=1.$$

The function around the origin is behaving very wildly.

I was thinking of breaking this guy up into a system of two first order ODE's and then solve, but I have no idea how to set this up. What method should I use to set up the system of ODE's?

If there is some other method rather than numerically solving a system of differential equations, please feel welcome to share. Thanks.

alt text

$\endgroup$
8
  • 2
    $\begingroup$ Is $f(1)=1$ really what you want? This will give division by zero on the right-hand side! Also, do you really have a boundary value problem (conditions at two different points $R=0$ and $R=1$), or did you mean an initial value problem with two conditions at one point (say $R=0$)? An IVP is relatively straightforward to integrate numerically, but a BVP is harder. $\endgroup$
    – Hans Lundmark
    Commented Oct 16, 2010 at 13:39
  • $\begingroup$ I am trying to replicate a graph of the solution from an old paper. I can show a picture of what the graph looks like but I am not sure how. Anyway the solution is for R between 0 and 1. At f(1)=1 at least from what the graphs shows. My other boundary condition is actually wrong (I will edit the main post), but f seems to be approaching positive infinity as R approached 0. $\endgroup$
    – theAngryPhysicist
    Commented Oct 16, 2010 at 21:43
  • $\begingroup$ Well I am not able to post the graph yet. I don't enough reputation points (i need 10). I might post a link later. $\endgroup$
    – theAngryPhysicist
    Commented Oct 16, 2010 at 21:54
  • $\begingroup$ More than posting the graph, where did this graph come from, and how did you construct your DE from said graph? (+1 so you can post the graph). $\endgroup$
    – J. M. ain't a mathematician
    Commented Oct 17, 2010 at 16:20
  • $\begingroup$ This is from an old research paper on Inertial Electrostatic Confinement, in English its a large vacuum tank with hot plasma inside. This equation describes the normalized potential of the plasma at different radii from the center of this spherical tank. $\endgroup$
    – theAngryPhysicist
    Commented Oct 19, 2010 at 15:29

2 Answers 2

Reset to default
1
$\begingroup$

The general method to reinterpret a higher-order ODE as a system of first order ODEs is to regard the derivatives of the unknown function as additional unknown functions. In your case, regard $f'$ as a new function $g$. Then the system of first order ODEs consists of two equations, the first being the original equation with $f'$ replaced by $g$ and $f''$ replaced by $g'$, and the second equation being $f'=g$.

$\endgroup$
4
  • $\begingroup$ Okay in that case I will end up with: g'(R)+(2/R)g(R)=(.7/R)((1/sqrt(h))-((0.3)/sqrt(1-h))), g=h' by letting f=h and f'=g. How will this effect the one boundary condition f(1)=1 I knew before the substitution? $\endgroup$
    – theAngryPhysicist
    Commented Oct 16, 2010 at 22:45
  • $\begingroup$ You don't need h. You just say f'=g and write the equation you have in terms of g',g, and f. It looks like yours with f instead of h. You still use your f(1)=1 (which divides by zero as stated before), but need a second boundary condition. $\endgroup$
    – Ross Millikan
    Commented Oct 17, 2010 at 0:52
  • $\begingroup$ Okay, could I use an approximation? $\endgroup$
    – theAngryPhysicist
    Commented Oct 19, 2010 at 15:39
  • $\begingroup$ For second-order ODEs you need two boundary conditions, otherwise you end up with a family of solutions (unless you're trying to pick a specific member using, say, optimization). Typically, you'll need a b.c. on f', which, after the transformation mentioned above, gives you a b.c. on g. To solve numerically the (stiff) ODE, the first thing to try is a Runge-Kutta method such as ode45 in Matlab. $\endgroup$
    – Dominique
    Commented Oct 25, 2011 at 22:51
0
$\begingroup$

I'm assuming you are interested in the region $R<1$; the region $R>1$ as shown in the graph.

Dominique is right, you do need a second constraint or boundary condition. More about that below.

Both the wicked behavior near 0 and the combination $f''+\frac{2}{R}f'$ argue for making the substitution $x = \frac{1}{R}$. This transforms the equation to $$ \frac{d^2f}{dx^2} = \frac{K_+}{x^3}\left( \sqrt{f} - \lambda_+\sqrt{1-f} \right) $$ where per the notation on the graph, $K_+$ is $0.7$ and $\lambda_+$ is $0.3$.

This equation is to be solved numerically for $x > 1$. In that region, the equation is quite well behaved, and Runge-Kutta will work well. Depending on how near to zero you need to get, you will have to integrate to some value $x_f > 1$. In order to get the solution down to $R=\epsilon$ you would integrate to $x_f = 1/\epsilon$.

And now the zinger: What is that second boundary condition?

If you are given $\frac{df}{dR}$ at $R=1$, this is an initial value problem starting from $x=1$, and that is straightforward. For example, the dashed curve on the graph looks to have $$\left. \frac{df}{dR} \right|_{R=1} = K_+ = 0.7$$

However, the main interesting feature on the graph is that $f(R=0.3) = 0$. Your graph contains multiple distinct dashed curves (which might be one reason you found it so hard to reproduce) since the dashed curve to the left of $0.3$ does not match the one to the right. This suggests that the boundary conditions used were $$ \left\{ \begin{array}{c} f(1) = 1 \\ f(\lambda_+) = 0 \end{array} \right.$$
Runge-Kutta is not ideal for a BVP but you can do an iterative solution by shooting for the fixed point at $x=1/0.3, f(x) = 0$, or there are relaxation methods that will work fine because the function is well behaved in this region.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .