在https://blog.csdn.net/weixin_42141390/article/details/110184743一文中，我们曾经讨论了欧拉法，龙格-库塔法也跟欧拉法一样，是用梯形的面积去替代积分的面积的一种方法。

欧拉法简介

设有微分方程：
dx(t)dt=f(x)\frac{dx(t)}{dt} = f(x) dtdx(t)=f(x)
x(t0)=x0x(t_0)=x_0x(t0)=x0已知。

我们对上述方程进行积分，积分区域为[t0,t1],Δt=t1−t0[t_0,t_1],\Delta t = t_1-t_0[t0,t1],Δt=t1−t0，可得：
x(t1)=x(t0)+∫t0t1f(t)dtx(t_1)=x(t_0)+\int_{t_0}^{t_1}f(t)dt x(t1)=x(t0)+∫t0t1f(t)dt
其中x(t1)x(t_1)x(t1)就是我们要求解的东西了。于是，就差积分怎么解决了。

对于欧拉法来说，积分是靠梯形面积来近似，如下所示：
∫t0t1f(t)dt=f(x(t0))(t1−t0)=f(x0)Δt\int_{t_0}^{t_1}f(t)dt = f(x(t_0))(t_1-t_0)=f(x_0)\Delta t ∫t0t1f(t)dt=f(x(t0))(t1−t0)=f(x0)Δt
带入上述方程即可得：
x(t1)=x(t0)+f(x0)Δtx(t_1)=x(t_0)+f(x_0)\Delta t x(t1)=x(t0)+f(x0)Δt
按照上式进行递推，即可得：
xk+1=xk+f(xk)Δtx_{k+1} = x_{k} + f(x_{k})\Delta t xk+1=xk+f(xk)Δt
其中 x0x_0x0 已知。问题就搞定了。

二阶龙格-库塔法

对于微分方程 dxdt=f(x)\frac{dx}{dt} = f(x)dtdx=f(x)
由于梯形的上底取得不合理，导致欧拉法存在一定的误差和不稳定性，于是需要对其进行改进。龙格库塔法的递推公式为：
xk=xk−1+Δt(λ1m1+λ1m2)x_{k} = x_{k-1} + \Delta t(\lambda_{1}m_1+\lambda_1 m_2) xk=xk−1+Δt(λ1m1+λ1m2)
其中：m1=f(xk−1)m_1 = f(x_{k-1})m1=f(xk−1)、m2=f(xk−1+m1⋅pΔt)m_2=f(x_{k-1}+m_1\cdot p\Delta t)m2=f(xk−1+m1⋅pΔt)，且 p⋅λ2=1/2,λ1+λ2=1p\cdot \lambda_2 = 1/2, \lambda_1+\lambda_2=1p⋅λ2=1/2,λ1+λ2=1

若：p=1,λ1=1/2,λ2=1/2p = 1, \lambda_1=1/2, \lambda_2=1/2p=1,λ1=1/2,λ2=1/2，则二阶龙格-库塔法等价于改进欧拉法。

高阶龙格-库塔法

我们可以对二阶龙格-库塔法进行推广，如三阶龙格-库塔法的递推公式为：
xk=xk−1+Δt6(m1+4m2+m3)x_{k} = x_{k-1} + \frac{\Delta t}{6}(m_1 + 4 m_2 + m_3) xk=xk−1+6Δt(m1+4m2+m3)
其中：m1=f(xk−1),m2=f(xk−1+1/2⋅m1Δt),m3=f(xk−1+(2m2−m1)Δt)m_1 = f(x_{k-1}), m_2 = f(x_{k-1}+1/2 \cdot m_1 \Delta t), m_3 = f(x_{k-1}+ (2m_2-m1)\Delta t)m1=f(xk−1),m2=f(xk−1+1/2⋅m1Δt),m3=f(xk−1+(2m2−m1)Δt)

也有四阶龙格库塔法（最常用）：
xk=xk−1+Δt6(m1+2m2+2m3+m4)x_{k} = x_{k-1} + \frac{\Delta t}{6}(m_1 + 2 m_2 + 2 m_3+m_4) xk=xk−1+6Δt(m1+2m2+2m3+m4)
其中：m1=f(xk−1),m2=f(xk−1+1/2⋅m1Δt),m3=f(xk−1+1/2⋅m2Δt),m4=f(xk−1+m3Δt)m_1 = f(x_{k-1}), m_2 = f(x_{k-1}+1/2 \cdot m_1 \Delta t), m_3 = f(x_{k-1}+ 1/2\cdot m_2 \Delta t), m_4 = f(x_{k-1}+m_3\Delta t)m1=f(xk−1),m2=f(xk−1+1/2⋅m1Δt),m3=f(xk−1+1/2⋅m2Δt),m4=f(xk−1+m3Δt)

就如我们在欧拉法这篇文章讲的那样，我们可以用泰勒公式来确定龙格库塔法的精确度，可以证明的是，多少阶的龙格库塔法，就有多少阶的精确度。

案例与代码

可见求解区域为 [0，1]，我们设置求解的步数为 100，也即 Δt=0.01s\Delta t = 0.01sΔt=0.01s，代码如下：

import numpy as np
import matplotlib.pyplot as plt
def f(x):return -20*x
def lgkt3(n,x0,t0,tn,f):x = [x0]t = [t0]for i in range(n):dt = (tn-t0)/ntk = t[i]+dtt.append(tk)m1 = f(x[i])*dtm2 = f(x[i]+m1*dt/2)m3 = f(x[i]+(2*m2-m1)*dt)xk = x[i]+dt/6*(m1+4*m2+m3)x.append(xk)return x,tdef lgkt4(n,x0,t0,tn,f):x = [x0]t = [t0]for i in range(n):dt = (tn-t0)/ntk = t[i]+dtt.append(tk)m1 = f(x[i])*dtm2 = f(x[i]+m1*dt/2)m3 = f(x[i]+m2*dt/2)m4 = f(x[i]+m3*dt)xk = x[i]+dt/6*(m1+2*m2+2*m3+m4)x.append(xk)return x,t
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = Falsefont1 = {'family' : 'SimHei',
'weight' : 'normal',
'size'   : 15,
}def myplot2(n,x0,t0,tn,f):x1,t = lgkt3(n,x0,t0,tn,f)x2,t = lgkt4(n,x0,t0,tn,f)plt.figure(figsize=(12,4))plt.subplot(1,2,1)plt.plot(t,x1,linewidth=3,color='r',label='3阶龙格库塔法')    plt.xlabel('t',fontsize=24)plt.ylabel('x',fontsize=24)plt.legend(prop=font1)plt.grid()plt.subplot(1,2,2)plt.plot(t,x2,linewidth=3,color='b',label='4阶龙格库塔法')plt.xlabel('t',fontsize=24)plt.ylabel('x',fontsize=24)plt.legend(prop=font1)plt.grid()if __name__ == '__main__':myplot2(10,1,0,1,f)

龙格库塔法求解微分方程相关推荐

龙格库塔法和欧拉法求解微分方程的比较
文章目录计算机如何理解连续系统的动态特性? 欧拉法求解微分方程龙格库塔法求解微分方程 MATLAB代码编写和仿真效果计算机如何理解连续系统的动态特性? 一般连续系统的动态特性可以由一个微分方程, ...
龙格库塔(runge-kutta，RK)法求解微分方程
求解微分方程的意思就是,已知导函数,求原函数. 先声明一点,欧拉法.中值法.龙哥库塔法求解微分方程,得出的结果不是表达式,而是一系列离散点. 一.欧拉法递推问题描述:已知y'=f(x,y),求y(x ...
四阶龙格库塔法的基本思想_利用龙格库塔法求解郎之万方程.doc
利用龙格库塔法求解郎之万方程.doc 利用龙格-库塔法求解朗之万方程1. 待解问题布朗颗粒是非常微小的宏观颗粒,其直径的典型大小为10-710-6m.颗粒不断受到液体介质分子的碰撞,在任一瞬间,一个颗 ...
Matlab求解微分方程数值解
有三种方法求解微分方程数值解: 欧拉法改进欧拉法龙格库塔法接下来用一个练习来对比这三种求解方法. 问题描述:用改进的Euler方法.MATLAB的ode45命令分别求下列初值问题的数值解,并画图 ...
matlab解二阶微分方程组,[微分方程组]急急急！用MATLAB按二阶龙格库塔法求解微分方程组，急用于毕业设计！...
急急急!用MATLAB按二阶龙格库塔法求解微分方程组,急用于毕业设计! 问题补充:今天才发现自己之前做的一点都不对,17号就交论文了,我傻了,急死了!求各位大侠帮帮忙.谢谢!要求解的微分方程如图所示. ...
AI攻破高数核心，1秒内精确求解微分方程、不定积分，性能远超Matlab
栗子鱼羊发自海边边量子位报道 | 公众号 QbitAI 大家都知道,AI (神经网络) 连加减法这样的简单算术都做不好: 可现在,AI已经懂得微积分,把魔爪伸向你最爱的高数了. 它不光会求 ...
龙格库塔法解微分方程组的matlab程序,MATLAB实例源码教程：龙格库塔法求解微分方程组源代码实例.doc...
MATLAB实例源码教程:龙格库塔法求解微分方程组源代码实例.doc MATLAB实例源码教程龙格库塔法求解微分方程组源代码实例题目用经典 Runge-Kutta方法求下列一阶微分方程组的近似解y1 ...
MATLAB从入门到精通-欧拉法与梯形法求解微分方程（含MATLAB源码）
前言以下是我为大家准备的几个精品专栏,喜欢的小伙伴可自行订阅,你的支持就是我不断更新的动力哟! MATLAB-30天带你从入门到精通 MATLAB深入理解高级教程(附源码) tableau可视化数据 ...
matlab龙格库塔法求通解,基于matlab及龙格库塔法求解布拉修斯方程.doc
基于matlab及龙格库塔法求解布拉修斯方程 Runge-Kutta法求解布拉修斯解摘要薄剪切层方程主要有三种解法,即相似解,非相似条件下对偏微分方程组的数值解和近似解.布拉修斯解是布拉修斯于19 ...

龙格库塔法求解微分方程

欧拉法简介

二阶龙格-库塔法

高阶龙格-库塔法

案例与代码

龙格库塔法求解微分方程相关推荐

最新文章

热门文章