Numerical solution of a (system of?) DE

Apr 2008
1,000
100
Teyateyaneng
I've got the 3 DE's:
(1)\(\displaystyle \ddot r= \frac{k}{m}(1-r)+r \dot \theta ^2 + r \dot \phi ^2 \sin (\theta) + g\cos (\theta)\).
(2)\(\displaystyle \ddot \phi = -2 \dot \phi \left [ \frac{\dot r}{r} - \dot \theta \cot (\theta) \right ]\).
(3)\(\displaystyle \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: \(\displaystyle r(0)=1\), \(\displaystyle \phi (0)=-\frac{\pi}{2}\), \(\displaystyle \theta (0)=\frac{\pi}{2}\).
And \(\displaystyle \dot r(0)=-\frac{1}{2}\), \(\displaystyle \dot \phi (0)=0\), \(\displaystyle \dot \theta (0)= \pi\).
One may assume that k, m and g are given. I want to get numerical values for \(\displaystyle r(t_i)\), \(\displaystyle \phi (t_i)\) and \(\displaystyle \theta (t_i)\) for a certain number of \(\displaystyle 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.
 

shawsend

MHF Hall of Honor
Aug 2008
903
379
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.
 
Apr 2008
1,000
100
Teyateyaneng
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 \(\displaystyle \omega=\dot \phi\), thus I'd have \(\displaystyle \dot w =-2 \omega \left [ \frac{\dot r}{r}- \dot \theta \cot (\theta) \right ]\). But this \(\displaystyle \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. >>
 
Last edited:

shawsend

MHF Hall of Honor
Aug 2008
903
379
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}]
 
  • Like
Reactions: arbolis