Help with using the Runge-Kutta 4th order method

  • #1
quasarLie
51
0

Homework Statement


I have this equations of motion, I have this equations of motion for a schwarchild black hole, I wish to use the 4th order Runge-Kutta method to solve them for a body falling to the black hole from a distance r0 and with L = 0. My problem is I am struggling to apply this method to my system of ODE's so that I can program a method that can solve any system of ODE's using the formulas for RK4, I would like for someone to please run through one step of the method, so I can understand it better Thank

Homework Equations


[/B]
\begin{eqnarray}
\frac{dt}{d\tau} &=& (1- \frac{1}{r})^{-1} \frac{E}{m c^{2}}\\
\frac{d\phi}{d\tau} &= & \dfrac{L}{mr^{2}} \\
(\dfrac{dr}{d\tau})^{2} &=& \frac{E^{2}}{m^{2} c^{4}}(1- \frac{R_{s}}{r}) \dfrac{L^{2}}{m^{2} c^{2} r^{2}}
\end{eqnarray}​

The Attempt at a Solution


Now I know that for a general 1st order ODE's the 4th order Runge-Kutta formula's:
\begin{equation}
y_{i+1}=y_i + \frac{1}{6}(k_1+2k_2+2k_3+k_4)
\end{equation}
with
\begin{equation}
k_1 = f(t_n,y_n) \\ k_2 = f(t_n+\frac{1}{2}h,y_n+\frac{h}{2}k_1) \\ k_3 = f(t_n+\frac{1}{2}h,y_n+\frac{h}{2}k_2) \\ k_4 = f(t_n+\frac{1}{2}h,y_n+\frac{h}{2}k_3)
\end{equation}
 
Last edited by a moderator:
Physics news on Phys.org
  • #2
You posted this originally in the Precalculus section -- this is hardly a precalc type of problem
quasarLie said:

Homework Statement


I have this equations of motion, I have this equations of motion for a schwarchild black hole, I wish to use the 4th order Runge-Kutta method to solve them for a body falling to the black hole from a distance r0 and with L = 0.
To the best of my knowledge, you can apply this method only to a single differential equation, not a system of differential equations.
quasarLie said:
My problem is I am struggling to apply this method to my system of ODE's so that I can program a method that can solve any system of ODE's using the formulas for RK4, I would like for someone to please run through one step of the method, so I can understand it better Thank
Which step are you uncertain of?
quasarLie said:

Homework Equations


[/B]
\begin{eqnarray}
\frac{dt}{d\tau} &=& (1- \frac{1}{r})^{-1} \frac{E}{m c^{2}}\\
\frac{d\phi}{d\tau} &= & \dfrac{L}{mr^{2}} \\
(\dfrac{dr}{d\tau})^{2} &=& \frac{E^{2}}{m^{2} c^{4}}(1- \frac{R_{s}}{r}) \dfrac{L^{2}}{m^{2} c^{2} r^{2}}
\end{eqnarray}​

The Attempt at a Solution


Now I know that for a general 1st order ODE's the 4th order Runge-Kutta formula's:
\begin{equation}
y_{i+1}=y_i + \frac{1}{6}(k_1+2k_2+2k_3+k_4)
\end{equation}
with
\begin{equation}
k_1 = f(t_n,y_n) \\ k_2 = f(t_n+\frac{1}{2}h,y_n+\frac{h}{2}k_1) \\ k_3 = f(t_n+\frac{1}{2}h,y_n+\frac{h}{2}k_2) \\ k_4 = f(t_n+\frac{1}{2}h,y_n+\frac{h}{2}k_3)
\end{equation}
 
  • #3
is this correct for the first equation
\begin{equation}
k_1=f(t_n,\tau_n) = (1- \frac{1}{r_n})^{-1} \frac{E}{m c^{2}}
\end{equation}
 
Last edited:
  • #4
Mark44 said:
To the best of my knowledge, you can apply this method only to a single differential equation, not a system of differential equations.

This is not correct. RK4 works fine for systems of ODEs. You just have to interpret the [itex]y_n[/itex] as multicomponent vectors. In your case, [itex]\tau[/itex] is your independent variable, and your [itex] y_n[/itex] s are:
[tex] y_n = (t_n(\tau), \phi_n(\tau), r_n(\tau))[/tex]
Here is a link which explains in more detail.
 
  • Like
Likes scottdave
  • #5
Mark44 said:
To the best of my knowledge, you can apply this method only to a single differential equation, not a system of differential equations.

phyzguy said:
This is not correct. RK4 works fine for systems of ODEs. You just have to interpret the yn as multicomponent vectors.
I had something like that in mind, but didn't say it. What I was trying to get across was that you couldn't just plug the whole system of equations into RK. You're solving each equation for yi separately.
 
  • #6
Mark44 said:
I had something like that in mind, but didn't say it. What I was trying to get across was that you couldn't just plug the whole system of equations into RK. You're solving each equation for yi separately.
Yes i want to solve each esuation for y
 
  • #7
From post #2:
Mark44 said:
Which step are you uncertain of?
 
  • #8
Mark44 said:
From post #2:
I want to know if this is correct for the second equation
\begin{equation}
k_{1}= \frac{L}{mr^{2}}\\
K_{2}= f(\phi(i)+\dfrac{h}{2},\dfrac{L}{mr_{i}^{2}}+ \frac{h}{2} \dfrac{L}{mr_{i}^{2}}\\
K_{3}= f(\phi_{i}+\frac{h}{2}, \dfrac{L}{mr_{i}^{2}}+\dfrac{h}{2}\dfrac{L}{mr_{i}^{2}})\\
K_{4}= f(\phi_{i}+h,\dfrac{L}{mr_{i}^{2}}+h\dfrac{L}{mr_{i}^{2}})
\end{equation}
 
  • #9
From post #1:
## \frac{dt}{d\tau} = (1- \frac{1}{r})^{-1} \frac{E}{m c^{2}} ##
## \frac{d\phi}{d\tau} = \frac{L}{mr^2} ##
## (\frac{dr}{d\tau})^{2} = \frac{E^{2}}{m^{2} c^{4}}(1- \frac{R_{s}}{r}) \frac{L^{2}}{m^{2} c^{2} r^{2}}##
I think you're going to have to do some work on the third equation before you can tackle the other two equations, since neither of the first two equations actually involves τ.
The third equation could be rewritten as ##\frac {dr}{d\tau} = \frac E {mc^2} \sqrt{1 - \frac {R_s}r} \frac L {mcr}##. Here I'm assuming that all of the squared quantities are nonnegative, so for instance, ##\sqrt{m^2} = m##.
If you solve this (manually), you will have r as a function of ##\tau##, and you can replace r in the first and second equations.
Then the first equation will have dt/dτ equal to an expression that actually involves τ, and similar for the second equation.
 
  • #10
Mark44 said:
From post #1:
## \frac{dt}{d\tau} = (1- \frac{1}{r})^{-1} \frac{E}{m c^{2}} ##
## \frac{d\phi}{d\tau} = \frac{L}{mr^2} ##
## (\frac{dr}{d\tau})^{2} = \frac{E^{2}}{m^{2} c^{4}}(1- \frac{R_{s}}{r}) \frac{L^{2}}{m^{2} c^{2} r^{2}}##
I think you're going to have to do some work on the third equation before you can tackle the other two equations, since neither of the first two equations actually involves τ.
The third equation could be rewritten as ##\frac {dr}{d\tau} = \frac E {mc^2} \sqrt{1 - \frac {R_s}r} \frac L {mcr}##. Here I'm assuming that all of the squared quantities are nonnegative, so for instance, ##\sqrt{m^2} = m##.
If you solve this (manually), you will have r as a function of ##\tau##, and you can replace r in the first and second equations.
Then the first equation will have dt/dτ equal to an expression that actually involves τ, and similar for the second equation.
Yes I have work to do with last equation to apply the RK4, i just want ti know if what i did is correct to use it in a C++ program to solve the quations. Thank you for your answers
 
  • #11
quasarLie said:
Yes I have work to do with last equation to apply the RK4, i just want ti know if what i did is correct to use it in a C++ program to solve the quations. Thank you for your answers
No, it doesn't look correct to me.
Here is your second equation: ## \frac{d\phi}{d\tau} = \frac{L}{mr^2} ##
To use Runge-Kutta, the equation needs to be in this form: ##\frac {d\phi}{d\tau} = f(\phi, \tau)##. That is, the right side of the equation has to be in terms of ##\phi## and ##\tau##. What you have in your second equation immediately above doesn't involve either of these variables. That's why I'm suggesting that you find a solution to the third equation, so that you can eventually get r in terms of ##\tau##, then substitute for r in the first two equations.
 
Last edited:
  • #12
Mark44 said:
No, it doesn't look correct to me.
Here is your second equation: ## \frac{d\phi}{d\tau} = \frac{L}{mr^2} ##
To use Runge-Kutta, the equation needs to be in this form: ##\frac {d\phi}{d\tau} = f(\phi, \tau)##. That is, the right side of the equation has to be in terms of ##\phi## and ##\tau##. What you have in your second equation immediately above doesn't involve either of these variables. That's why I'm suggesting that you find a solution to the third equation, so that you can eventually get r in terms of ##\tau##, then substitute for r in the first two equations.
This is what I did for the RK2. Do you think that it's false ? (it's a python program)
tRK2[i+1] = t + dtau * (t + dtau/2 + (1/(1-(Rs/R)) * E) + (dtau/2)*(1/(1-(Rs/R)) * E))
phiRK2[i+1] = phi + dtau * (phi + dtau/2 + (L/(m*R**2)) + (dtau/2)*(L/(m*R**2)))
RRK2[i+1] = R + dtau * (R + dtau/2 + V + (dtau/2)*V)
VRK2[i+1] = VRK2 - dtau * (VRK2 + dtau/2 + ((Rs/(2*R**2)) + (3/2)*((Rs*L**2 - R*L**2)/(R**4))) + (dtau/2)*((Rs/(2*R**2)) + (3/2)*((Rs*L**2 - R*L**2)/(R**4))))
 
  • #13
quasarLie said:
This is what I did for the RK2. Do you think that it's false ? (it's a python program)
Python:
 tRK2[i+1] = t[I] + dtau * (t[I] + dtau/2 + (1/(1-(Rs/R[I])) * E) + (dtau/2)*(1/(1-(Rs/R[I])) * E))
        phiRK2[i+1] = phi[I] + dtau * (phi[I] + dtau/2 + (L/(m*R[I]**2)) + (dtau/2)*(L/(m*R[I]**2)))
        RRK2[i+1] = R[I] + dtau * (R[I] + dtau/2 + V[I] + (dtau/2)*V[I])
        VRK2[i+1] = VRK2[I] - dtau * (VRK2[I] + dtau/2 + ((Rs/(2*R[I]**2)) + (3/2)*((Rs*L**2 - R[I]*L**2)/(R[I]**4))) + (dtau/2)*((Rs/(2*R[I]**2)) + (3/2)*((Rs*L**2 - R[I]*L**2)/(R[i]**4))))[/i][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/I][/QUOTE]
Your program is unreadable -- all expressions such as t get converted to BBCode for italics and disappear. Also, when you post code, you should use code tags -- see https://www.physicsforums.com/threa...-for-nearly-200-programming-languages.774303/ for more information. Using code tags retains the indentation, which is critical for Python code.
Second, I've invested a fair amount of time in perusing what you've done for RK4. Why don't you proceed with that, rather than pursue a different direction?
 
  • #14
Mark44 said:
Your program is unreadable -- all expressions such as t get converted to BBCode for italics and disappear. Also, when you post code, you should use code tags -- see https://www.physicsforums.com/threa...-for-nearly-200-programming-languages.774303/ for more information. Using code tags retains the indentation, which is critical for Python code.
Second, I've invested a fair amount of time in perusing what you've done for RK4. Why don't you proceed with that, rather than pursue a different direction?

Thank you for the information. and I just want to show you what i did for RK2 (which is correct) because I want to do the same thing with RK4 this is why I wrote the K_i
 
  • #15
quasarLie said:
Yes I have work to do with last equation to apply the RK4, i just want ti know if what i did is correct to use it in a C++ program to solve the quations. Thank you for your answers

Besides what Mark44 said in post #9, another possibiity is that
$$ \frac {dr}{d\tau} = -\frac E {mc^2} \sqrt{1 - \frac {R_s}r} \frac L {mcr}$$
That will give different dynamics from the choice
$$ \frac {dr}{d\tau} = \frac E {mc^2} \sqrt{1 - \frac {R_s}r} \frac L {mcr},$$
so you need to try to figure out which is the relevant form.
 

Similar threads

Back
Top