Doing 3 second order differential equations numerically

  • Thread starter WraithM
  • Start date
  • #1
32
0
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 eachother. :/

These are the equations:

[tex]\ddot{t} = \frac{-2M \dot{r} \dot{t}}{r(r-2M)}[/tex]

[tex]\ddot{\phi} = \frac{-2\dot{r}\dot{\phi}}{r}[/tex]

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

Dots mean differentiation with respect to [tex]\tau[/tex].

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

Any help would be greatly appreciated.
 

Answers and Replies

  • #2
32
0
Actually, this could be in C++ too. C or C++ works.
 
  • #3
240
2
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
[tex]\vec{x}=(x_{i}), i=1,2,3,4,5,6[/tex],

and define a six-component vector-field

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

So your equations looks like this:

[tex]\frac{d}{d\tau}\vec{x}=\vec{F}(\vec{x})[/tex]

The components of the state variable will be:

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

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:
[tex]F_{1}=x_{4}[/tex]
[tex]F_{2}=x_{5}[/tex]
[tex]F_{3}=x_{6}[/tex]
[tex]F_{4}=\frac{-2M x_{6} x_{4}}{x_{3} (x_{3}-2M)}[/tex]
[tex]F_{5}=\frac{-2x_{6} x_{5}}{x_{3}}[/tex]
[tex]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}[/tex]

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

[tex]\dot{\vec{x}}=\vec{F}(\vec{x})[/tex]

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.
 
  • #4
70
0
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.
 
  • #5
32
0
Thank you elibj123! Your method worked perfectly. I replicated my results in Mathematica, and my results in C almost perfectly fit the Mathematica curve.
 

Related Threads on Doing 3 second order differential equations numerically

  • Last Post
Replies
3
Views
2K
  • Last Post
Replies
2
Views
2K
  • Last Post
Replies
5
Views
2K
Replies
1
Views
2K
Replies
6
Views
3K
Replies
1
Views
5K
  • Last Post
Replies
9
Views
941
Replies
2
Views
8K
Replies
5
Views
1K
Replies
6
Views
3K
Top