Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Doing 3 second order differential equations numerically

  1. Feb 11, 2010 #1
    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.
     
  2. jcsd
  3. Feb 12, 2010 #2
    Actually, this could be in C++ too. C or C++ works.
     
  4. Feb 12, 2010 #3
    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.
     
  5. Feb 12, 2010 #4
    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.
     
  6. Feb 13, 2010 #5
    Thank you elibj123! Your method worked perfectly. I replicated my results in Mathematica, and my results in C almost perfectly fit the Mathematica curve.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Doing 3 second order differential equations numerically
Loading...