在https://blog.csdn.net/weixin_/article/details/一文中，我们曾经讨论了欧拉法，龙格-库塔法也跟欧拉法一样，是用梯形的面积去替代积分的面积的一种方法。
欧拉法简介
我们对上述方程进行积分，积分区域为 [ t 0 , t 1 ] , Δ t = t 1 − t 0 [t_0,t_1],\Delta t = t_1-t_0 [t0,t1],Δt=t1−t0，可得：
x ( t 1 ) = x ( t 0 ) + ∫ t 0 t 1 f ( t ) d t x(t_1)=x(t_0)+\int_{t_0}^{t_1}f(t)dt x(t1)=x(t0)+∫t0t1f(t)dt
其中 x ( t 1 ) x(t_1) x(t1)就是我们要求解的东西了。于是，就差积分怎么解决了。
二阶龙格-库塔法
对于微分方程 d x d t = f ( x ) \frac{dx}{dt} = f(x) dtdx=f(x)
由于梯形的上底取得不合理，导致欧拉法存在一定的误差和不稳定性，于是需要对其进行改进。龙格库塔法的递推公式为：
x k = x k − 1 + Δ t ( λ 1 m 1 + λ 1 m 2 ) x_{k} = x_{k-1} + \Delta t(\lambda_{1}m_1+\lambda_1 m_2) xk=xk−1+Δt(λ1m1+λ1m2)
其中： m 1 = f ( x k − 1 ) m_1 = f(x_{k-1}) m1=f(xk−1)、 m 2 = f ( x k − 1 + m 1 ⋅ 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 = 1 p\cdot \lambda_2 = 1/2, \lambda_1+\lambda_2=1 p⋅λ2=1/2,λ1+λ2=1
若： p = 1 , λ 1 = 1 / 2 , λ 2 = 1 / 2 p = 1, \lambda_1=1/2, \lambda_2=1/2 p=1,λ1=1/2,λ2=1/2，则二阶龙格-库塔法等价于改进欧拉法。
高阶龙格-库塔法
就如我们在欧拉法这篇文章讲的那样，我们可以用泰勒公式来确定龙格库塔法的精确度，可以证明的是，多少阶的龙格库塔法，就有多少阶的精确度。
案例与代码
可见求解区域为 [0，1]，我们设置求解的步数为 100，也即 Δ t = 0.01 s \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)/n tk = t[i]+dt t.append(tk) m1 = f(x[i])*dt m2 = 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,t def lgkt4(n,x0,t0,tn,f): x = [x0] t = [t0] for i in range(n): dt = (tn-t0)/n tk = t[i]+dt t.append(tk) m1 = f(x[i])*dt m2 = 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'] = False font1 = {
'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) 
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/220732.html原文链接：https://javaforall.net
