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

  • 2
    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.2010-10-16
  • 0
    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.2010-10-16
  • 0
    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.2010-10-16
  • 0
    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).2010-10-17
  • 0
    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.2010-10-19
  • 0
    Since you talk about spherical tank, I guess you have the homogeneous Neumann boundary condition at 0. It would help if you post the original 3D equation that must hold inside the tank.2010-12-24
  • 0
    You might find [this blog post](http://blogs.mathworks.com/loren/2013/06/10/from-symbolic-differential-equations-to-their-numeric-solution/) from The MathWorks helpful.2014-12-01
  • 0
    possible duplicate of [How to reduce higher order linear ODE to a system of first order ODE?](http://math.stackexchange.com/questions/501745/how-to-reduce-higher-order-linear-ode-to-a-system-of-first-order-ode) – note that there is no difference between linear and nonlinear in the answer given. This question has been answered many times on this site so I'm sure you can find more examples if you search.2014-12-01

2 Answers 2

1

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

  • 0
    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?2010-10-16
  • 0
    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.2010-10-17
  • 0
    Okay, could I use an approximation?2010-10-19
  • 0
    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.2011-10-25
0

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.