四阶龙格-库塔法——以2022年高教杯数学建模A题为例

一切都是瞬息

于 2024-07-17 13:58:50 发布

阅读量1.3k

点赞数 16

文章标签：数学建模 python

本文链接：https://blog.csdn.net/add70093ms/article/details/140472404

版权

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

文章目录

前言
一、原理
二、示例

前言

在微分方程的数值求解当中，四阶龙格－库塔法（Runge-Kutta methods）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。本文以2022年高教杯数学建模A题为例，详细讲解四阶龙格－库塔法的原理及其python实现。

一、原理

假设我们求解以下微分方程
$\left\{\begin{aligned}&\frac{dy}{dt}=f(y,t) \\&y(0)=y_0\end{aligned}\right.$
由中值定理可得
$y_{n+1}=y_{n}+f(y_{\varepsilon },t_{\varepsilon })\Delta t$
其中， $f(y_{\varepsilon },t_{\varepsilon })$ 为区间 $t_n,t_{n+1})$ 上的平均斜率。四阶龙格－库塔法的基本思想，就是找一个容易计算且与 $f(y_{\varepsilon },t_{\varepsilon })$ 相差极小的值来代替 $f(y_{\varepsilon },t_{\varepsilon })$ 进行运算。

四阶龙格－库塔法的表达式为
$y_{n+1}=y_{n}+\frac{\Delta t}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})$
其中
$\begin{aligned} &k_1=f(y_n,t_n) \\ &k_2=f\left(y_n+k_1\frac {\Delta t}2,t_n+\frac {\Delta t}2\right) \\ &k_3=f\left(y_n+k_2\frac {\Delta t}2,t_n+\frac {\Delta t}2\right) \\ &k_4=f\left(y_n+k_3\Delta t,t_n+\Delta t\right) \end{aligned}$
四阶龙格－库塔法的总体截断误差： $\propto\Delta t^4$

二、示例

2022年高教杯数学建模A题求解问题1.(1)的微分方程为（推导过程略）
$\begin{aligned} (m_{f}+m_{a})\frac{\mathrm{d}^{2}z_{f}}{\mathrm{d}t^{2}}& =f\cos\omega t-c_{a}\frac{\mathrm{d}z_{f}}{\mathrm{d}t}-k_{f}z_{f}-c_{p}\left(\frac{\mathrm{d}z_{f}}{\mathrm{d}t}-\frac{\mathrm{d}z_{v}}{\mathrm{d}t}\right)-k_{s}(z_{f}-z_{v}), \\ m_v\frac{\mathrm{d}^2z_v}{\mathrm{d}t^2}& =-c_p\left(\frac{\mathrm{d}z_v}{\mathrm{d}t}-\frac{\mathrm{d}z_f}{\mathrm{d}t}\right)-k_s(z_v-z_f), \end{aligned}$
初始条件为
$t=0:\quad z_f(0)=z_v(0)=\frac{\mathrm{d}z_f}{\mathrm{d}t}(0)=\frac{\mathrm{d}z_v}{\mathrm{d}t}(0)=0.$
将其化为一阶常微分方程
$\begin{aligned} &\frac{\mathrm{d}z_1}{\mathrm{d}t} =v_{1}, \\ &\frac{\mathrm{d}v_1}{\mathrm{d}t} =\frac1{m_f+m_a}[f\cos\omega t-c_av_1-k_fz_1-c_p(v_1-v_2)-k_s(z_1-z_2)], \\ &\frac{\mathrm{d}z_2}{\mathrm{d}t} =v_2, \\ &\frac{\mathrm{d}v_2}{\mathrm{d}t} =-\frac1{m_v}[c_p(v_2-v_1)+k_s(z_2-z_1)]. \end{aligned}$
根据四阶龙格－库塔法，我们可以列出 $z_1$ 的递推公式
$z1_{n+1}=z1_{n}+\frac{\Delta t}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})$
其中
$\begin{aligned} &f(v1_v)=v1_n\\ &g\left(z1_n,v1_n,z2_n,v2_n,t_n\right)=\frac1{m_f+m_a}[f\cos\omega t-c_av1_n-k_fz1_n-c_p(v1_n-v2_n)-k_s(z1_n-z2_n)] \end{aligned}$
$\begin{aligned} &k_1=f\left(v1_n\right) \\ &k_2=f\left(v1_n+l_1\frac {\Delta t}2\right) \\ &k_3=f\left(v1_n+l_2\frac {\Delta t}2\right)\\ &k_4=f\left(v1_n+l_3\Delta t\right) \end{aligned}$

$l_1,l_2,l_3$ 是 $v 1$ 的增长率，由 $g\left(z1,v1,z2,v2,t\right)$ 得到。同理，我们可以得到其它三个变量的递推公式。

四阶龙格－库塔法写出各变量的递推式后，我们可以开始编辑python代码。以下是一个示例

import numpy as np
import matplotlib.pyplot as pl
pl.rcParams['font.sans-serif'] = ['SimHei']
pl.rcParams['axes.unicode_minus'] = False
# 相关参数
omega = 1.4005
mf = 4866
ma = 1335.535
mv = 2433
f = 6250
ca = 656.3616
kf = 1025 * 9.8 * np.pi * 1 ** 2
cp = 10000
ks = 80000

# 步长和迭代次数
dt = 0.2
n = 600

# 初始条件
z1 = 0
z2 = 0
v1 = 0
v2 = 0
t = 0

z1_list = [0]
z2_list = [0]
v1_list = [0]
v2_list = [0]
t_list = [0]

# v1的增长率
def function1(z1,v1,z2,v2,t):
    dv1 = 1/(mf+ma)*(f*np.cos(omega*t)-ca*v1-kf*z1-cp*(v1-v2)-ks*(z1-z2))
    return dv1

# v2的增长率
def function2(z1,v1,z2,v2,t):
    dv2 = -1/mv*(cp*(v2-v1)+ks*(z2-z1))
    return dv2

# 四阶龙格-库塔法
for i in range(n):
    dz1_1 = v1
    dz2_1 = v2
    dv1_1 = function1(z1,v1,z2,v2,t)
    dv2_1 = function2(z1,v1,z2,v2,t)

    dz1_2 = v1 + dv1_1 * dt / 2
    dz2_2 = v2 + dv2_1 * dt / 2
    dv1_2 = function1(z1 + dz1_1 * dt / 2, v1 + dv1_1 * dt / 2, z2 + dz2_1 * dt / 2, v2 + dv2_1 * dt / 2, t + dt / 2)
    dv2_2 = function2(z1 + dz1_1 * dt / 2, v1 + dv1_1 * dt / 2, z2 + dz2_1 * dt / 2, v2 + dv2_1 * dt / 2, t + dt / 2)

    dz1_3 = v1 + dv1_2 * dt / 2
    dz2_3 = v2 + dv2_2 * dt / 2
    dv1_3 = function1(z1 + dz1_2 * dt / 2, v1 + dv1_2 * dt / 2, z2 + dz2_2 * dt / 2, v2 + dv2_2 * dt / 2, t + dt / 2)
    dv2_3 = function2(z1 + dz1_2 * dt / 2, v1 + dv1_2 * dt / 2, z2 + dz2_2 * dt / 2, v2 + dv2_2 * dt / 2, t + dt / 2)

    dz1_4 = v1 + dv1_3 * dt
    dz2_4 = v2 + dv2_3 * dt
    dv1_4 = function1(z1 + dz1_3 * dt, v1 + dv1_3 * dt, z2 + dz2_3 * dt, v2 + dv2_3 * dt, t + dt)
    dv2_4 = function2(z1 + dz1_3 * dt, v1 + dv1_3 * dt, z2 + dz2_3 * dt, v2 + dv2_3 * dt, t + dt)

    z1 = z1 + dt / 6 * (dz1_1 + 2 * dz1_2 + 2 * dz1_3 + dz1_4)
    z2 = z2 + dt / 6 * (dz2_1 + 2 * dz2_2 + 2 * dz2_3 + dz2_4)
    v1 = v1 + dt / 6 * (dv1_1 + 2 * dv1_2 + 2 * dv1_3 + dv1_4)
    v2 = v2 + dt / 6 * (dv2_1 + 2 * dv2_2 + 2 * dv2_3 + dv2_4)
    t = t + dt

    z1_list.append(z1)
    z2_list.append(z2)
    v1_list.append(v1)
    v2_list.append(v2)
    t_list.append(t)

print(t_list[50],z1_list[50],v1_list[50],z2_list[50],v2_list[50])
print(t_list[100],z1_list[100],v1_list[100],z2_list[100],v2_list[100])
print(t_list[200],z1_list[200],v1_list[200],z2_list[200],v2_list[200])
print(t_list[300],z1_list[300],v1_list[300],z2_list[300],v2_list[300])
print(t_list[500],z1_list[500],v1_list[500],z2_list[500],v2_list[500])

fig = pl.figure(figsize=(10,4))
ax1 =fig.add_subplot(1,2,1)
ax2 =fig.add_subplot(1,2,2)

ax1.plot(t_list, z1_list, 'b-', label='浮子',linewidth=1.0)
ax2.plot(t_list, z2_list, 'r-', label='振子',linewidth=1.0)

pl.subplots_adjust(left=0.15,bottom=0.1,top=0.9,right=0.95, \
                   hspace=0.35,wspace=0.3)
ax1.set_ylabel('z/m')
ax1.set_xlabel('T/s')
ax1.set_xlim(0,120)
ax1.set_ylim(-1.0,1.0)
ax2.set_ylabel('z/m', )
ax2.set_xlabel('T/s')
ax2.set_xlim(0,120)
ax2.set_ylim(-1.0,1.0)
ax1.legend(loc='upper right')
ax2.legend(loc='upper right')
pl.show()