1、龙格-库塔法简介
龙格—库塔法是一种在工程上应用广泛的高精度单步算法,其中包括著名的欧拉法,用于数值求解微分方程。 由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。 在各种龙格—库塔法当中有一个方法十分常用,该方法主要是在已知方程导数和初值信息,利用计算机仿真时应用,省去求解微分方程的复杂过程,该方法被称为“RK4”或者就是“龙格—库塔法”。
2、龙格-库塔法基本思想
在了解龙格—库塔法之前先回顾一下欧拉法。
2.1 欧拉法基本原理
若已知上述方程组某点坐标,以及该点的斜率值,则可以使用欧拉法对 y(x) 进行迭代的近似求解。求解分析:已知 x = x0 时 y = y0 y=y ,且斜率为 y′( x0) = f ( x0 , y0 ) ,那么就可以用过该点且斜率为 y′( x0) 的直线近似求解 x0 附近一点 x1 的函数值。
构建直线:
那么,在 x = x1 时,求得:
可以求得 y = y ( x1 ),重复以上步骤,可以求得 y(x2) , y(x3) , ⋯ , y (xn) 。如果等间隔取 x 值,即 x1−x0= x2− x1 = x3 − x2 = ⋯ = xn − xn−1 = h,也即 xn = x0 + n ⋅ h 可以得到
总结为:
该公式就是欧拉公式。
2.2 龙格-库塔基本原理
龙格-库塔算法是对欧拉法的改进,都是使用过 ( xk , yk) 的直线来近似下一点的函数值 yk+1 ,区别在于斜率的选择。在欧拉法中,直接使用该点的切线 y′(xk) 当做直线的斜率,这样会造成对于 y(xk+1) 的求解一定存在误差,除非 y(x) 为直线时,结果才是准确的。基于上面的考虑,需要解决的便是如何提高 y(xk+1) 求解的精度;或者说在知道 ( xk , yk) 的情况下,如何确定直线斜率确保恰好经过待估测点 y(xk+1) 。答案是利用拉格朗日中值定理。
拉格朗日中值定理如下:
即在 ( xk , xk+1) 区间内存在一点 ϵ (xk< ϵ < xk+1) ,使得等式 f (xk+1) − f (xk) = f ′( ϵ ) (xk+1 − xk) 成立。
所以准确的直线方程为:
如果等间距取值,即 xk+1 − xk = h ,则有
其中,θ ∈ ( 0 , 1 ) ,而当 θ = 0 时就是欧拉法。
.
现在,原理已经明确了,但如何确定 y′(ϵ) 的取值?这就涉及到算法的精度问题。考虑到直接使用端点处 y′( xk ) 的斜率不准确,那么取两端点斜率的平均准确度肯定就提升了一点,即y′(ϵ) 取
为了方便书写,将 y′(xk) = κ1, y′(xk+1) = κ2 ,即用 κ 表示斜率。则上式这可改写为:
根据已知条件方程 y′(x) = f(x,y) ,同理可得:
其中 f(x,y) 是已知的,而 f(xk+1 , yk+1) 中的 yk+1 是未知的,若用欧拉法来估计此处的 yk+1,即有:
当然这是不准确的,若用两点的平均斜率来提高精,则有:
值得注意的是,对于两端点斜率取平均是出于习惯,但实际应用中是对两斜率进行加权,取平均只是习惯性的把两者权重取相等而已,因此更一般描述为:
.
根据上述斜率取平均值的思路,进行递推取 n 个点,获得更多点加权来估算 f′(ϵ) 的值所得结果精度肯定更高。现在,假设在 (xk, xk+1) 区间取 n 个点,对应点为 xk+α1h , xk+α2h , xk+α3h , ⋯ , xk+αnh , 且 α ∈ [ 0 , 1 ] ,在每个点的斜率记为 κ1, κ2, ⋯ , κn ,对应加权为 λ1 , λ2, ⋯ , λn , 那么有
则有:
λ 为自取值,对于 κ 有:
通常 α1 取0,即对应 xk 点;此时,κ1 = f (xk, yk) 。同样,对于其余式子 f 自变量中的 y(xk+ αnh ) 使用欧拉法来求,即有:
总结以上,即得到龙格-库塔基本表达式:
其中:λ1, λ2, ⋯ , λn 为对应采样点斜率权值,满足 λ1 + λ2 + ⋯ + λn = 1 ;
α1, α2, ⋯ , αn 为采样点位置比率,满足 0 < α1 < α2 < ⋯ < αn
.
3、 案例:四阶龙格-库塔法求解微分方程
本次案例使用的微分方程如下:
本次案例使用的是常用的 经典R-K方法,如下:
文章来源:https://www.toymoban.com/news/detail-616006.html
.
运行代码文章来源地址https://www.toymoban.com/news/detail-616006.html
# lgkt
import numpy as np
def RK(x0, y0, h, n, func):
"""
x0: 初始点x坐标
y0: 初始点y坐标
h: 步长
n: 迭代次数
func: 事先定义好的f(x,y)
"""
x = np.linspace(x0, x0 + (n - 1) * h, num=n)
y = np.zeros_like(x)
y[0] = y0
for i in range(n - 1):
k1 = func(x[i], y[i])
k2 = func(x[i] + h / 2, y[i] + h * k1 / 2)
k3 = func(x[i] + h / 2, y[i] + h * k2 / 2)
k4 = func(x[i] + h, y[i] + h * k3)
y[i + 1] = y[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
return x, y
if __name__ == '__main__':
# y'=f(x,y)=2x/(3y^2)
f = lambda x, y: 2 * x / (3 * y * y)
X, Y = RK(0, 1, h=0.4, n=5, func=f)
print("X=", X)
print("Y=", Y)
到了这里,关于龙格-库塔方法学习笔记的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!