对于常微分方程,RK方法速度快,精度高,代码简单,是最为实用的数值方法之一。RK方法很简单,类似梯形法,RK法也是根据前一步点的值推算后一步点。具体算法见以下链接
https://wenku.baidu.com/view/3d7e77450a4c2e3f5727a5e9856a561252d32184.html?from=search&smallflow20190502=1
例:求解如下微分方程
精确解)
显式方法中步长不能太大,否则结果不可靠。在代码中28对步长进行了判断。
例子中步长要小于0.066。 读者可取n=5, 10,100对比下。
函数RK计算四个系数,fxy是方程右端的表达式。
import math
import numpy as np
from numpy import *
import matplotlib.pyplot as plt
import time
def RK(y0, a, b, n):#输入y0,x的区间【a,b】以及等分数
h = (b-a)/n
y = np.zeros(n)
y[0] = y0
for i in range(1, n, 1): #从1到n
x0 = a+(i-1)*h # 这里对应上一步的x0
k1 = fxy(x0, y0, h)
k2 = fxy(x0+h/2., y0+h/2.*k1, h)
k3 = fxy(x0+h/2., y0+h/2.*k2, h)
k4 = fxy(x0+h, y0+h*k3, h)
y0 = y0+h/6.*(k1+2.*k2+2.*k3+k4)
y[i] = y0
i = i+1
return y
def fxy(x, y, h): #被积函数写在这里
lanmb=-x*x*y #lanmb为正数的时候不用判断
f = lanmb*y #这里需要判断步长是否收敛。表达式dy/dx=lanmb*y
if (lanmb*h < -2):
print('h should smaller than ', abs(2/lanmb), h)# 收敛判断条件
return f
start = time.clock()
a0 = 0.
b0 = 1.5
y0 = 3. #y的初始值
n = 80
yy = RK(y0, a0, b0+ (b0-a0)/n, n+1)# 这里是闭区间,开区间不需要加 (b0-a0)/n
xx = np.arange(a0, b0+ (b0-a0)/n, (b0-a0)/n)
print(xx[-1],yy[-1])#打印最后一个点
yyy=3./(1+xx**3) #精确解 测试用
plt.figure(1) #画图
print(xx.shape, yy.shape)
plt.plot(xx, yy)
plt.scatter(xx,yyy) #精确解
delta=np.sum(abs(yyy-yy))
print(delta)
end = time.clock() #看一下所用时间
print('time=',end-start)
plt.show()
最后结果:散点是精确解,曲线是数值解