# Help with using the Runge-Kutta 4th order method

## 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:

Mark44
Mentor
You posted this originally in the Precalculus section -- this is hardly a precalc type of problem

## 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}

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:
phyzguy
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 $y_n$ as multicomponent vectors. In your case, $\tau$ is your independent variable, and your $y_n$ s are:
$$y_n = (t_n(\tau), \phi_n(\tau), r_n(\tau))$$
Here is a link which explains in more detail.

• scottdave
Mark44
Mentor
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 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.

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

Mark44
Mentor
From post #2:
Which step are you uncertain of?

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}

Mark44
Mentor
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.

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

Mark44
Mentor
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:
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))))

Mark44
Mentor
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[i] 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 persue a different direction?

Your program is unreadable -- all expressions such as t[i] 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 persue 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

Ray Vickson
Homework Helper
Dearly Missed
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.