python龙格库塔_龙格库塔积分法

Runge--Kutta积分法是由Euler方法改进得来

As with the forward-difference derivative, the error in Euler’s rule is O(h2), which is equivalent to ignoring the effect of acceleration on the position of a projectile, but including it on the projectile’s velocity.

在许多计算中，Euler方法的精度(具有1阶精度)是不够的。这个时候，我们就需要一个更好的算法。

栗子

原例中的matlab代码

改造后的Python code

import numpy as np

import matplotlib.pyplot as plt

import matplotlib

matplotlib.rcParams['font.family']='STSong'

fig=plt.figure(figsize=(12,8), dpi=100)

N=20

h=0.25

X,Y=[0],[2]

def f(x,y):

return -x*y**2

y_n=2

for i in range(N):

x_n=i*h

k_1=f(x_n,y_n)

k_2=f(x_n+0.5*h,y_n+0.5*h*k_1)

k_3=f(x_n+0.5*h,y_n+0.5*h*k_2)

k_4=f(x_n+h,y_n+h*k_3)

y_n+=1/6*h*(k_1+2*k_2+2*k_3+k_4)

X.append(x_n+h)

Y.append(y_n)

plt.plot(X,Y,'r:')

x_n,y_n

0.25 1.8823080259996157

0.5 1.5998962064862998

0.75 1.2799477682362133

1.0 1.0000271050738332

1.25 0.7805556108821389

1.5 0.6154593936835502

1.75 0.4923742156717609

2.0 0.4000542886135871

2.25 0.3299396220734262

2.5 0.2758952271913286

2.75 0.23360233153202317

3.0 0.20001998162311824

3.25 0.17298862395359263

3.5 0.15095575873910502

3.75 0.13278993601441305

4.0 0.11765498312076311

4.25 0.10492446251452127

4.5 0.0941229086206177

4.75 0.08488497716443959

5.0 0.07692668513821321

# rk4Algor.py: algorithm for Delta y, input y, f; do NOT modify

import numpy as np

def rk4Algor(t, h, N, y, f):

k1=np.zeros(N); k2=np.zeros(N); k3=np.zeros(N); k4=np.zeros(N)

k1 = h*f(t,y)

k2 = h*f(t+h/2.,y+k1/2.)

k3 = h*f(t+h/2.,y+k2/2.)

k4 = h*f(t+h,y+k3)

y=y+(k1+2*(k2+k3)+k4)/6.

return y

y" = -100y-2y'+ 10 sin(3t)

# rk4Call.py: 4th-O Runge-Kutta that calls rk4Algor

# Here for ODE y" = -100y-2y'+ 10 sin(3t)

from rk4Algor import rk4Algor

from numpy import *

import numpy as np, matplotlib.pyplot as plt

Tstart = 0.; Tend = 10.; Nsteps = 100 # Initialization

tt=[]; yy=[]; yv=[]; y = zeros((2), float) #

y[0] = 3.; y[1] = -5. # Initial position & velocity

t = Tstart; h = (Tend-Tstart)/Nsteps;

def f(t, y): # Force (RHS) function

fvector = zeros((2), float)

fvector[0] = y[1]

fvector[1] = -100.*y[0]-2.*y[1] + 10.*sin(3.*t)

return fvector

while (t < Tend):

tt.append(t) # Time loop

if ((t + h) > Tend): h = Tend - t # Last step

y = rk4Algor(t, h, 2, y, f) #

yy.append(y[0])

yv.append(y[1])

t = t + h

fig=plt.figure()

plt.subplot(111)

plt.plot(tt,yy,'r')

plt.title('Position versus')

plt.xlabel('t')

plt.ylabel('y')

fig1=plt.figure()

plt.subplot(111)

plt.plot(tt,yv)

plt.title('Velocity versus time')

plt.xlabel('t')

plt.ylabel('y')

plt.show