RK4 二阶微分方程 PYTHON

Question

我有一个用 RK4 求解的方程：

任何人都可以提供强大的解决方案吗？

我有一个带有参数的代码：

k = 33
m = 0.05
K_imp = 100
x_d = 0.05
x_0 = 0.1
v_0 = 0.0
h = 0.01
ddx_0=0
t_end = 10
update_interval=10
t = np.arange(0, t_end, h)
xpoints = []
upoints = []
step = 0

#Equation to solve: 
def F(x,u,t):
    return -x
for i in t:
    xpoints.append(x)
    upoints.append(u)
    
    m1=h*u
    k1 = h*F(y, u, x)  #(x, v, t)
    
    m2 = h*(u + 0.5*k1)
    k2 = h*F(y+0.5*m1, u+0.5*k1, x+0.5*h)

    m3 = h*(u + 0.5*k2)
    k3 = h*F(y+0.5*m2, u+0.5*k2, x+0.5*h)

    m4 = h*(u + k3)
    k4 = h*F(y+m3, u+k3, x+h)

    y += (m1 + 2*m2 + 2*m3 + m4)/6
    u += (k1 + 2*k2 + 2*k3 + k4)/6
    step=step+1
if step % update_interval == 0:
           x_jm1 = x[i]
           step=0

我无法编写代码，我需要一些帮助，x_d 是一段时间后所需的解决方案。

Answer 1

@vitya，你的方程的RHS没有意义，所以我用x替换了x_j-1。那么至少你可以写下确切的解决方案。

除非您在方程的 LHS 中添加阻尼项，否则解永远不会趋向于最终值 x_d。

下面是基本的龙格-库塔解决方案。

import numpy as np
import math


def rk4( x, y, dx, f ):
    dy1 = dx * f( x           , y             )
    dy2 = dx * f( x + 0.5 * dx, y + 0.5 * dy1 )
    dy3 = dx * f( x + 0.5 * dx, y + 0.5 * dy2 )
    dy4 = dx * f( x +       dx, y +       dy3 )
    return x + dx, y + ( dy1 + 2.0 * dy2 + 2.0 * dy3 + dy4 ) / 6.0


# Equation parameters
k = 33
m = 0.05
K_imp = 100
x_d = 0.05
x_0 = 0.1
v_0 = 0.0
def deriv( t, x ):
    return np.array( [ x[1], ( K_imp * ( x_d - x[0] ) - k * x[0] ) / m ] )


# Exact solution
w = math.sqrt( ( k + K_imp ) / m )
A = x_0 - x_d * K_imp / ( k + K_imp )
B = v_0 / w
print( "Period = {:.4f}".format( 2.0 * math.pi / w ) )
def exact( t ):
    return A * math.cos( w * t ) + B * math.sin( w * t ) + x_d * K_imp / ( k + K_imp )


# Computational parameters
dt = 0.002         # timestep
nt = 200           # number of timesteps
ofreq = 10         # output frequency


# Initialise
t, x = 0.0, np.array( [ x_0, v_0 ] )
fmt = " {:11.6f}"

# Successive timesteps
for i in range( nt + 1 ):
    if i % ofreq == 0: print( (3*fmt).format( t, x[0], exact( t ) ) )
    t, x = rk4( t, x, dt, deriv )

列是 t、x、exact_x（可能）：

Period = 0.1218
    0.000000    0.100000    0.100000
    0.020000    0.069641    0.069641
    0.040000    0.008103    0.008102
    0.060000   -0.024743   -0.024743
    0.080000    0.003062    0.003062
    0.100000    0.064464    0.064464
    0.120000    0.099723    0.099724
    0.140000    0.074534    0.074534
    0.160000    0.013405    0.013404
    0.180000   -0.024190   -0.024191
    0.200000   -0.001673   -0.001673
    0.220000    0.059049    0.059050
    0.240000    0.098896    0.098897
    0.260000    0.079100    0.079100
    0.280000    0.018921    0.018920
    0.300000   -0.023091   -0.023091
    0.320000   -0.006060   -0.006059
    0.340000    0.053444    0.053445
    0.360000    0.097526    0.097527
    0.380000    0.083298    0.083298
    0.400000    0.024603    0.024602