Doing 3 second order differential equations numerically

[WraithM]

I want to write a program in C that's going to solve 3 second order differential equations, and I have little to no experience in solving differential equations in C (however, I know C well enough). I'm having a bit of trouble.

First of all, I don't really know what method I'm going to use to solve these differential equations. I know how to do Euler's method, and Runge-Kutta and those things, but perhaps there's a math library out there for C that works magic with differential equations. I don't really care how these equations get solved. I just care that they get solved, but I'm more than willing to write some code myself to solve them. If anybody has any suggestions, that'd be extremely appreciated.

Also, I don't really know how I'm going to solve 3 equations at once. If I write my own differential equation solver, I'm going to substitute and make them first order, then solve again, but I don't really know how to deal with equations that depend on each other. :/

These are the equations:

\ddot{t} = \frac{-2M \dot{r} \dot{t}}{r(r-2M)}

\ddot{\phi} = \frac{-2\dot{r}\dot{\phi}}{r}

\ddot{r} = \frac{M\dot{r}^2}{r(r-2M)} - \frac{M(r-2M)}{r^3}\dot{t}^2 + (r - 2M) \dot{\phi}^2

Dots mean differentiation with respect to \tau.

How do you solve 3 differential equations that all depend upon each other at once (numerically)?

Any help would be greatly appreciated.

 

[WraithM]

Actually, this could be in C++ too. C or C++ works.

 

[elibj123]

I don't know about C++ but software which was design for mathematic analysis (especially Matlab) surely will have a lot of ODE/PDE tools.

In regards to the algorithm itself, it's not much different from one ODE.
And firstly it will be more simple to bring it down to a set of six first order ODE rather than three second-order.
What you do is define six component state variable
\vec{x}=(x_{i}), i=1,2,3,4,5,6,

and define a six-component vector-field

\vec{F}=(F_{i}(\vec{x})), i=1,2,3,4,5,6

So your equations looks like this:

\frac{d}{d\tau}\vec{x}=\vec{F}(\vec{x})

The components of the state variable will be:

x_{1}=t; x_{2}=\phi; x_{3}=r; x_{4}=\dot{t}; x_{5}=\dot{\phi}; x_{6}=\dot{r};

The first components of F will contain the information that the derivative of the first three variables, are actually the last three variables. The last three component of F will relate the derivatives of x4,5,6 to the equations themselves:
F_{1}=x_{4}
F_{2}=x_{5}
F_{3}=x_{6}
F_{4}=\frac{-2M x_{6} x_{4}}{x_{3} (x_{3}-2M)}
F_{5}=\frac{-2x_{6} x_{5}}{x_{3}}
F_{6}=\frac{M x^{2}_{6}}{x_{3} (x_{3}-2M)}-\frac{M ( x_{3}-2M)}{x^{3}_{3}} x^{2}_{4} + (x_{3}-2M)x^{2}_{5}

So we've reduced the equations to the forms of:

\dot{\vec{x}}=\vec{F}(\vec{x})

Now forget the vector signs, and just treat it like a first order ODE, which you know how to solve using Euler's method or Runge-Kutta.

For example in C\++, define a struct\class that will represent a six-component vector, with some basic arithmetic operations.

 

[kosovtsov]

Your ODE system can be solved exactly. Please use the following Maple code

ode1:=diff(t(tau),tau,tau)=-(2*M*diff(r(tau),tau)*diff(t(tau),tau))/(r(tau)*(r(tau)-2*M));

ode2:=diff(phi(tau),tau,tau)=-(2*diff(r(tau),tau)*diff(phi(tau),tau))/(r(tau));

ode3:=diff(r(tau),tau,tau)=(M*diff(r(tau),tau)^2)/(r(tau)*(r(tau)-2*M))-(M*(r(tau)-2*M))/r(tau)^3*diff(t(tau),tau)^2+(r(tau)-2*M)*diff(phi(tau),tau)^2;
ans:=[dsolve({ode1,ode2,ode3})];

The answer is lengthy enough to be displayed here.

 

[WraithM]

Thank you elibj123! Your method worked perfectly. I replicated my results in Mathematica, and my results in C almost perfectly fit the Mathematica curve.