Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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),$$


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

share|cite|improve this question
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. – Hans Lundmark Oct 16 '10 at 13:39
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. – Shagster_84 Oct 16 '10 at 21:43
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. – Shagster_84 Oct 16 '10 at 21:54
+1 to give you the reputation for posting the graph. – whuber Oct 17 '10 at 15:04
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). – J. M. Oct 17 '10 at 16:20

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$.

share|cite|improve this answer
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? – Shagster_84 Oct 16 '10 at 22:45
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. – Ross Millikan Oct 17 '10 at 0:52
Okay, could I use an approximation? – Shagster_84 Oct 19 '10 at 15:39
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. – Dominique Oct 25 '11 at 22:51

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.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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