使用ODEINT与Runge-Kutta-4th

Question

我正在建模一个耦合的弹簧质量系统：两个带有M1和M2的对象，并与弹簧常数K1和K2耦合，并带有阻尼D1和D2。

Method1：使用odeint进行

cookbook-script

来求解微分方程。

# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

def vectorfield(w, t, p):
    """
    Defines the differential equations for the coupled spring-mass system.

    Arguments:
        w :  vector of the state variables:
                  w = [x1,y1,x2,y2]
        t :  time
        p :  vector of the parameters:
                  p = [m1,m2,k1,k2,b1,b2]
    """
    x1, y1, x2, y2 = w
    m1, m2, k1, k2, b1, b2 = p

    f = [y1,
         (-d1 * y1 - k1 * x1  + k2 * (x2 - x1 )+ d2 * (y2-y1)) / m1,
         y2,
         (-d2 * (y2-y1) - k2 * (x2 - x1 )) / m2]

    return f

# Parameter values
tau=1.02 
f_1=0.16 
f_2=tau*f_1

m1 = 2000000 # mass1
m2 = 20000 # mass2

# Spring constants
k1 = m1*pow(2*np.pi*f_1,2)
k2 = m2*pow((2*np.pi*f_2),2)


# damping

d1 = (0.04/2/np.pi)*2*pow(k1*m1,0.5)



d_d1=6000
l_p=9.81/pow(2*np.pi*f_2,2)
b=0.3
d2=d_d1*pow((l_p-b)/l_p,2)

# Initial conditions
# x1 and x2 are the initial displacements; y1 and y2 are the initial velocities
x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 100.0
numpoints = 5001 

# Create the time samples for the output of the ODE solver.
# I use a large number of points, only because I want to make
# a plot of the solution that looks nice.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

# Pack up the parameters and initial conditions:
p = [m1, m2, k1, k2, d1, d2]
w0 = [x1, y1, x2, y2]

# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
              atol=abserr, rtol=relerr)

#    convert zip object to tuple object
temp=tuple(zip(t, wsol[:,0],wsol[:,1],wsol[:,2],wsol[:,3]))
#    convert tulpe object to list then to dataframe
df=pd.DataFrame(list(map(list,temp)))

# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1') 


legend = ax.legend(loc='upper right', shadow=None, fontsize='small')
ax.set_xlabel('time [s]', fontdict=None, labelpad=None, loc='center')
ax.set_ylabel('pos [m] or V [m/s]', fontdict=None, labelpad=None, loc='center')

Method2：使用Runge-Kutta-4th Order（由我自己编程）来求解微分方程

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.pyplot as plt




# Parameter values
tau=1.02 
f_1=0.16 
f_2=tau*f_1

m1 = 2000000 # mass1
m2 = 20000 # mass2

# Spring constants
k1 = m1*pow(2*np.pi*f_1,2)
k2 = m2*pow((2*np.pi*f_2),2)


# damping
d1 = (0.04/2/np.pi)*2*pow(k1*m1,0.5)



d_d1=6000
l_p=9.81/pow(2*np.pi*f_2,2)
b=0.3
d2=d_d1*pow((l_p-b)/l_p,2)

def system1(x1,y1,x2,y2):
    return (-d1 * y1 - k1 * x1 + k2 * (x2 - x1 )+ d2 * (y2-y1)) / m1 

def system2(x1, y1, x2, y2):
    return (-d2 * (y2-y1) - k2 * (x2 - x1)) / m2


def runge_kutta_4_sys1(f, x1, y1, x2, y2, h):
    k1 = f(x1,y1,x2,y2)
    k2 = f(x1,y1+k1*h/2,x2,y2)
    k3 = f(x1,y1+k2*h/2,x2,y2)
    k4 = f(x1,y1+k3*h,x2,y2)
 
    

    temp1= y1
    temp2= y1+k1*h/2
    temp3= y1+k2*h/2
    temp4= y1+k3*h
    
    
    displace_solution = x1 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6
    velocity_solution = y1 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6
    acc_solution = f(x1,y1,x2,y2)
 
    return displace_solution, velocity_solution, acc_solution

def runge_kutta_4_sys2(f, x1, y1, x2, y2, h):
    k1 = f(x1,y1,x2,y2)
    k2 = f(x1,y1,x2,y2+k1*h/2)
    k3 = f(x1,y1,x2,y2+k2*h/2)
    k4 = f(x1,y1,x2,y2+k3*h)
 
    

    temp1= y2
    temp2= y2+k1*h/2
    temp3= y2+k2*h/2
    temp4= y2+k3*h
    
    
    displace_solution = x2 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6
    velocity_solution = y2 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6
    acc_solution = f(x1,y1,x2,y2)
 
    return displace_solution, velocity_solution, acc_solution

x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0


h = 0.02
numpoints = 5000
time = 0

temp=0

df = pd.DataFrame(index=range(1+numpoints),columns=range(7))

df.iloc[0]=[time,x1,y1,temp,x2,y2,temp]

for i in range(numpoints):
    x1, y1, acc1 = runge_kutta_4_sys1(system1, x1, y1, x2, y2, h)
    x2, y2, acc2 = runge_kutta_4_sys2(system2, x1, y1, x2, y2, h)
    time=time+h
    df.iloc[i+1]=[time,x1,y1,acc1,x2,y2,acc2]
    df.iloc[i,3]=acc1
    df.iloc[i,6]=acc2

# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1') 


legend = ax.legend(loc='upper right', shadow=None, fontsize='small')
ax.set_xlabel('time [s]', fontdict=None, labelpad=None, loc='center')
ax.set_ylabel('pos [m] or V [m/s]', fontdict=None, labelpad=None, loc='center')    

方法1和方法之间的解决方案显然有偏差。

i我假设偏差是由于不准确的runge-kutta融合所致。我的假设正确吗？如何设置或提高runge-kutta-methode的准确性？

我尝试在Internet中搜索以使用Runger-Kutta方法找到可能的限制和准确性，而无需成功。

是在评论中首先观察到的@lastchance，您的RK4不正确，因为它不会同时将方法应用于所有四个方程。我也不了解RWO Velocity_solution方程背后的基本原理。

将其RK4方法应用于四个方程的系统。出于教学原因，我明确地添加了它的“光荣”复杂性。对于稍微复杂的系统，将不得不添加抽象，以允许在循环中计算K因子组。

def system1(x1, y1, x2, y2):
    return (-d1 * y1 - k1 * x1 + k2 * (x2 - x1) + d2 * (y2 - y1)) / m1

def system2(x1, y1, x2, y2):
    return (-d2 * (y2 - y1) - k2 * (x2 - x1)) / m2

def runge_kutta_4(f1, f2, x1, y1, x2, y2, h):
    k1x1 = y1
    k1x2 = y2
    k1y1 = f1(x1, y1, x2, y2)
    k1y2 = f2(x1, y1, x2, y2)

    k2x1 = y1 + h * k1y1 / 2
    k2x2 = y2 + h * k1y2 / 2
    k2y1 = f1(x1 + h * k1x1 / 2, y1 + h * k1y1 / 2, x2 + h * k1x2 / 2, y2 + h * k1y2 / 2)
    k2y2 = f2(x1 + h * k1x1 / 2, y1 + h * k1y1 / 2, x2 + h * k1x2 / 2, y2 + h * k1y2 / 2)

    k3x1 = y1 + h * k2y1 / 2
    k3x2 = y2 + h * k2y2 / 2
    k3y1 = f1(x1 + h * k2x1 / 2, y1 + h * k2y1 / 2, x2 + h * k2x2 / 2, y2 + h * k2y2 / 2)
    k3y2 = f2(x1 + h * k2x1 / 2, y1 + h * k2y1 / 2, x2 + h * k2x2 / 2, y2 + h * k2y2 / 2)

    k4x1 = y1 + h * k3y1
    k4x2 = y2 + h * k3y2
    k4y1 = f1(x1 + h * k3x1, y1 + h * k3y1, x2 + h * k3x2, y2 + h * k3y2)
    k4y2 = f2(x1 + h * k3x1, y1 + h * k3y1, x2 + h * k3x2, y2 + h * k3y2)

    x1_next = x1 + h * (k1x1 + 2 * k2x1 + 2 * k3x1 + k4x1) / 6
    x2_next = x2 + h * (k1x2 + 2 * k2x2 + 2 * k3x2 + k4x2) / 6
    y1_next = y1 + h * (k1y1 + 2 * k2y1 + 2 * k3y1 + k4y1) / 6
    y2_next = y2 + h * (k1y2 + 2 * k2y2 + 2 * k3y2 + k4y2) / 6

    acc1_next = f1(x1_next, y1_next, x2_next, y2_next)
    acc2_next = f2(x1_next, y1_next, x2_next, y2_next)

    return x1_next, x2_next, y1_next, y2_next, acc1_next, acc2_next

x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

h = 0.02
numpoints = 5000
time = 0

temp1 = system1(x1, y1, x2, y2)
temp2 = system1(x1, y1, x2, y2)

df = pd.DataFrame(index=range(1 + numpoints), columns=range(7))

df.iloc[0] = [time, x1, y1, temp1, x2, y2, temp2]

for i in range(numpoints):
    x1, x2, y1, y2, acc1, acc2 = runge_kutta_4(system1, system2, x1, y1, x2, y2, h)
    time = time + h
    df.iloc[i + 1] = [time, x1, y1, acc1, x2, y2, acc2]
    df.iloc[i, 3] = acc1
    df.iloc[i, 6] = acc2

Answer 1