Solving 2nd Order DEs with 4th Order RK Method

In summary, there are two methods discussed for rewriting equations in the second order case and implementing the fourth order Runge-Kutta method. However, there may be issues with the second method as it may not give the correct result for test cases and may not be the best choice for calculating a trajectory for an object.
  • #1
Arman777
Insights Author
Gold Member
2,168
192
Homework Statement
Create a trajectory for the comets orbit
Relevant Equations
##\ddot{x} = -GMx/r^3## and ##\ddot{y} = -GMy/r^3## for ##r = \sqrt{x^2 + y^2}##
In second order case we should rewrite the equation in terms of 2 first order DE's. So I wrote,

$$dx/dt = wx$$ $$dwx/dt = -GMx/r^3$$ and $$dy/dt = wy$$, $$dwy/dt = -GMy/r^3$$

Now I guess there's two ways to do it in 4th order RK method. I would either do it component by component or just in once. However when I do it component by component I couldn't uptade the distance .

Version1:
from numpy import arange, array
from math import sqrt
from pylab import plot, show, xlabel, ylabel

G = 6.67408 * 10 ** -11 #m^3kg^^-1s^-2
M = 1.989 * 10 ** 30 # kg
t_i = 0
t_f = 100
N = 10000
h = (t_f - t_i) / N

def f(r, t):
    x = r[0] # dx/dt = wx
    wx = r[1] #dwx/dt = -GMx/r^3
    y = r[2] # dx/dt = wx
    wy = r[3] #dwx/dt = -GMx/r^3
    f_x = wx
    f_wx = (- G * M * x) / sqrt(x**2 + y**2)**3
    f_y = wy
    f_wy = (- G * M * y) / sqrt(x**2 + y**2)**3
    return array([f_x, f_wx, f_y, f_wy], float)

t_data = arange(t_i, t_f, h)
x_points = []
y_points = []

r = array([4*10**12, 0, 0, 500], float) # [x, vx, y, vy]for t in t_data:
    x_points.append(r[0])
    y_points.append(r[2])
    k1 = h * f(r, t)
    k2 = h * f(r + 0.5 * k1, t + 0.5 * h)
    k3 = h * f(r + 0.5 * k2, t + 0.5 * h)
    k4 = h * f(r + k3, t + h)
    r += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)plot(x_points, y_points)
show()

or

Version2:
from numpy import arange, array
from pylab import plot, show, xlabel, ylabel
from math import sqrt

G = 6.67408 * 10 ** -11 #m^3kg^^-1s^-2
M = 1.989 * 10 ** 30 # kg
t_i = 0
t_f = 2000
N = 10000
h = (t_i - t_f) / N

def fx(rx, ry, t):
    x = rx[0] # dx/dt = wx
    wx = rx[1] #dwx/dt = -GMx/r^3
    y = ry[0]
    f_x = wx
    f_wx = (- G * M * x) / sqrt(x**2 + y**2) **3
    return array([f_x, f_wx], float)

def fy(rx, ry, t):
    y = ry[0] # dx/dt = wx
    wy = ry[1] #dwx/dt = -GMx/r^3
    x = rx[0]
    f_y = wy
    f_wy = (- G * M * y) / sqrt(x**2 + y**2) **3
    return array([f_y, f_wy], float)

t_data = arange(t_i, t_f, N)
x_points = []
y_points = []

rx = array([4*10**12, 0], float) #position and velocity respectively
ry = array([0, 500], float)

# for the x-component
for t in t_data:
    x_points.append(rx[0])
    y_points.append(ry[0])
    k1 = h * fx(rx, ry, t)
    k2 = h * fx(rx + 0.5 * k1, ry, t + 0.5 * h)
    k3 = h * fx(rx + 0.5 * k2, ry, t + 0.5 * h)
    k4 = h * fx(rx + k3, ry, t + h)
    k1 = h * fy(rx, ry, t)
    k2 = h * fy(rx, ry + 0.5 * k1, t + 0.5 * h)
    k3 = h * fy(rx, ry + 0.5 * k2, t + 0.5 * h)
    k4 = h * fy(rx, ry+ k3, t + h)
    rx += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)
    ry += 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)

plot(x_points, y_points)
show()

I am open to suggestions
 
Last edited:
Physics news on Phys.org
  • #2
I'm not sure what you mean by "couldnt uptade the distance".
One problem I see in the second method: You update rx and then calculate the change in y using the old ry but the new rx.
 
  • #3
mfb said:
I'm not sure what you mean by "couldnt uptade the distance".
One problem I see in the second method: You update rx and then calculate the change in y using the old ry but the new rx.
So I Should change the position of the rx. So you think that my code is correct ?

I changed my version 2 code.
 
  • #4
Does it give the correct result for test cases?
 
  • #5
mfb said:
Does it give the correct result for test cases?
Not really, I guess it also does not matter at this point. It seems that the RK method is not the best choice to calculate a trajectory for an object.
 
  • Skeptical
Likes BvU

What is a 2nd Order Differential Equation?

A 2nd order differential equation is a mathematical equation that involves a function, its first and second derivatives, and independent variable. It can be written in the form of y'' = f(x,y,y').

What is the 4th Order Runge-Kutta Method?

The 4th Order Runge-Kutta (RK4) method is a numerical method used to solve ordinary differential equations. It is an iterative method that uses a set of equations to approximate the solution of a differential equation at a given point. It is a popular method due to its high accuracy and efficiency.

How do you solve a 2nd Order DE with the 4th Order RK Method?

To solve a 2nd order differential equation using the 4th Order RK Method, you need to first convert it into a system of two 1st order differential equations. Then, you can use the RK4 method to iteratively find the solution at each time step. This involves calculating four approximations of the solution at each step and taking their weighted average.

What are the advantages of using the 4th Order RK Method?

The 4th Order RK Method has several advantages over other numerical methods for solving differential equations. It is highly accurate and can handle stiff equations. It also has a variable step size, which allows for more efficient computation. Additionally, it is a widely used and well-understood method in the field of numerical analysis.

What are some common applications of solving 2nd Order DEs with the 4th Order RK Method?

The 4th Order RK Method is commonly used in various fields, such as physics, engineering, and economics, to solve differential equations that model real-world phenomena. Some specific applications include modeling population growth, motion of objects under the influence of forces, and electrical circuits.

Similar threads

  • Engineering and Comp Sci Homework Help
Replies
6
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
3
Views
1K
  • Programming and Computer Science
Replies
15
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
3
Views
1K
  • Engineering and Comp Sci Homework Help
Replies
1
Views
2K
  • Engineering and Comp Sci Homework Help
Replies
4
Views
3K
  • Engineering and Comp Sci Homework Help
Replies
4
Views
3K
  • Advanced Physics Homework Help
Replies
1
Views
929
  • Atomic and Condensed Matter
Replies
3
Views
766
  • Engineering and Comp Sci Homework Help
Replies
7
Views
2K
Back
Top