Adams隐式4阶方法解常微分方程，python实现

Adams隐式4阶方法解常微分方程，由4阶Runge-Kutta方法提供初值，隐式方法比显式复杂一些，主要是因为需要解方程。这里使用弦截法解微分方程。

import math
import numpy as np  
import matplotlib.pyplot as plt
import time
def Secant(y3,x3,y2,x2,y1,x1,y0,x0,h):                  #这里y3 表示未知数
    eps=1.e-12
    i = 0                                      # 
    dx0 = 0.1                  # dx0用于差商，这里是初始值
    while(abs(dx0)>eps and i < 20):  #控制循环次数
        dfx0 = (Fun(y3+dx0,x3,y2,x2,y1,x1,y0,x0,h)-Fun(y3,x3,y2,x2,y1,x1,y0,x0,h))/dx0  #x点处的差商
        if dfx0==0:
            print('dfx0=0,y3=',y3)
            return y3
        y31 = y3-Fun(y3,x3,y2,x2,y1,x1,y0,x0,h)/dfx0             #计算新的x
        dx0 = y31-y3              #新的dx，用于求导 ，这行和牛顿法略有不同
        y3 = y31
        i = i+1
    if abs(Fun(y3,x3,y2,x2,y1,x1,y0,x0,h))>10e-4:             #判断f(x)是不是根，如果不是返回99999,然后主程序里面可以将这个值过滤掉
        return 999999        
    print(y3) 
    return y3

def RK(y0,  a, b, n):#   RK 法，计算前几个值输入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 Fun(y3,x3,y2,x2,y1,x1,y0,x0,h):                       #这里是要解的方程，dy/dx=x*x*y*y，需要输入前三步的计算结果
    fn3=-x3*x3*y3*y3
    fn2=-x2*x2*y2*y2
    fn1=-x1*x1*y1*y1
    fn0=-x0*x0*y0*y0
    f =  y3-y2-h/24.*(9.*fn3+19.*fn2-5.*fn1+fn0)
    return f

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

def Adams(a0,b0,y00,h):     #Adams 法接着计算后面的值，这里采用四阶Adams 隐式公式
    x = np.arange(a0, b0+h, h) #闭区间，所以用b00+h
    y = np.zeros(x.size)
    y[0] = y00 
    n = 3
    y0= RK(y[0],  a0, a0+3.*h, n)    #Rk法计算前三个值
    y[0:3] = y0
    for i in range(n, x.size,1):
        y[i] = Secant(y[i-1],x[i],y[i-1],x[i-1],y[i-2],x[i-2],y[i-3],x[i-3],h) # y[n-1]作为迭代初值
    return y 

start = time.clock()     # 计算时间
a0 = 0.                  # 定义域和步长
b0 = 1.5
y0 = 3.  #初始条件
h = 0.1
xx = np.arange(a0, b0+h, h)

yy = Adams(a0,b0,y0,h)     # 调用Adams方法

yyy = 3./(1+xx**3)         #解析解
delta = np.sum(abs(yyy-yy))  #计算误差
print(delta)

end = time.clock()
print('time=',end-start)

plt.figure(1)
plt.plot(xx, yy)
plt.scatter(xx,yyy)
plt.show()

最后通过比较误差发现实际上Runge-Kutta显式方法的误差更小一些。这个例子里面隐式方法更稳定的优势没体现出来