# Numerical solution of a (system of?) DE

• May 26th 2010, 05:49 PM
arbolis
Numerical solution of a (system of?) DE
I've got the 3 DE's:
(1) $\ddot r= \frac{k}{m}(1-r)+r \dot \theta ^2 + r \dot \phi ^2 \sin (\theta) + g\cos (\theta)$.
(2) $\ddot \phi = -2 \dot \phi \left [ \frac{\dot r}{r} - \dot \theta \cot (\theta) \right ]$.
(3) $\ddot \theta =\dot \phi ^2 \sin (\theta) \cos (\theta) -\frac{2 \dot r \dot \theta}{r} - \frac{g}{r} \sin (\theta)$.
With the initial conditions: $r(0)=1$, $\phi (0)=-\frac{\pi}{2}$, $\theta (0)=\frac{\pi}{2}$.
And $\dot r(0)=-\frac{1}{2}$, $\dot \phi (0)=0$, $\dot \theta (0)= \pi$.
One may assume that k, m and g are given. I want to get numerical values for $r(t_i)$, $\phi (t_i)$ and $\theta (t_i)$ for a certain number of $t_i$.
My problem is that they're second order DE's, while Euler's method and Runge-Kutta's work for first order DE. So I think I should reduce the order of the equations I have. How should I tackle this problem? A numerical integration? If so, I guess the error will be terribly big because I'd start to use Euler's method from an already "erroneous" function.
So I'm stuck at solving these equation. Any help is greatly appreciated as usual.
• May 27th 2010, 11:42 AM
shawsend
I see nothing wrong with just plugging those into Mathematica's NDSolve and going as long as you specify the constants.

You know how to do that right? Just plug them in as is, exactly as written and unless you run into singular points during the integration, there should be no problems.
• May 27th 2010, 02:09 PM
arbolis
Quote:

Originally Posted by shawsend
I see nothing wrong with just plugging those into Mathematica's NDSolve and going as long as you specify the constants.

You know how to do that right? Just plug them in as is, exactly as written and unless you run into singular points during the integration, there should be no problems.

Hey thank you for helping. But isn't NDSolve useful for first order DE's, like mentioned in the help file of Mathematica? I think it's a good idea if Mathematica can solve the problem for me, but as I need to use Euler and then a Runge-Kutta method, I'd use Mathematica only to confirm my result I get with Fortran.
Edit: I've found a similar problem on the internet but less complicated in which they reduced the order of the DE this way: by setting $\omega=\dot \phi$, thus I'd have $\dot w =-2 \omega \left [ \frac{\dot r}{r}- \dot \theta \cot (\theta) \right ]$. But this $\dot r$ really bothers me.
Edit2: Hmm, I've tried
Code:

s = NDSolve[{phi''[t] + 2*phi'*(r'/r - theta'*Cot[theta]) == 0,   phi[0] == -Pi/2, y'[0] == 0}, phi, {t, 0, 30}]
, but I get
Code:

NDSolve::dvnoarg: The function phi^\[Prime] appears with no arguments. >>
• May 27th 2010, 03:15 PM
shawsend
You can solve any order ODE, PDE that is reasonably solvable in Mathematica: Notice I specify a function as r[t] and theta[t] and so on. Remember to cut/paste it then change it to Standard Form via Cell/Convert To/Standard Form so it's easier to read. This is however Mathematica 7. If you have an earlier version there may be downward incompatibilities.

Code:

k = 1; m = 1; g = 1; eqn1 = Derivative[2][r][t] ==     (k/m)*(1 - r[t]) +     r[t]*Derivative[1][\[Theta]][t]^2*       Sin[\[Theta][t]] + g*Cos[\[Theta][t]]; eqn2 = Derivative[2][\[Phi]][t] ==     -2*Derivative[1][\[Phi]][t]*     (Derivative[1][r][t]/r[t] -       Derivative[1][\[Theta]][t]*Cot[\[Theta][t]]); eqn3 = Derivative[2][\[Theta]][t] ==     Derivative[1][\[Phi]][t]^2*Sin[\[Theta][t]]*       Cos[\[Theta][t]] -     (2*Derivative[1][r][t]*       Derivative[1][\[Theta]][t])/r[t] -     (g/r[t])*Sin[\[Theta][t]]; mysol = NDSolve[{eqn1, eqn2, eqn3,     r[0] == 1, Derivative[1][r][0] ==     -2^(-1), \[Phi][0] == -Pi/2,     Derivative[1][\[Phi]][0] == 0,     \[Theta][0] == Pi/2, Derivative[1][\[Theta]][       0] == Pi}, {r, \[Phi], \[Theta]}, {t, 0, 1}] Plot[Evaluate[{r[t], \[Phi][t], \[Theta][t]} /.     mysol], {t, 0, 1}, PlotStyle ->   {Red, Blue, Green}]