4阶显式Runge－Kutta法解常微分方程的通用程序--python实现

对于常微分方程，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()

最后结果：散点是精确解，曲线是数值解