常微分方程初值问题的数值解法

调包

import math
import numpy as np
import matplotlib.pyplot as plt

显示欧拉法

fStr为函数str名

#显式欧拉法
def EulerExplicit(x0,y0,h,fStr):
    xn = x0
    yn = y0
    n = 0
    ns = [n]
    xs = [xn]
    ys = ['%.8f'%yn]
    while n < 50:
        n += 1
        xn += h
        yn = yn + globals()[fStr](xn,yn) * h
        ns.append(n)
        xs.append('%.2f'%xn)
        ys.append('%.8f'%yn)
    return (ns,xs,ys)

ns为迭代次数，xs为x值列表，ys为y值列表

隐式欧拉法

#隐式欧拉法
def EulerImplicit(x0,y0,h,fStr):
    xn = x0
    yn = y0
    n = 0
    ns = [n]
    xs = [xn]
    ys = ['%.8f'%yn]
    while n < 50:
        n += 1
        xn +=h
        yp = yn + globals()[fStr](xn,yn) * h
        yn = yn + globals()[fStr](xn,yp) * h
        ns.append(n)
        xs.append('%.2f'%xn)
        ys.append('%.8f'%yn)
    return (ns,xs,ys)

欧拉改进法

#欧拉改进法
def EulerImproved(x0,y0,rlimX,h,fStr):
    xn = x0
    yn = y0
    n = 0
    ns = [n]
    xs = [xn]
    ys = ['%.8f'%yn]
    while True:
        yp = yn + globals()[fStr](xn,yn) * h
        xn += h
        n += 1
        if xn > rlimX:
            return (ns,xs,ys)
        yc = yn + globals()[fStr](xn,yp) * h
        yn = (yp + yc) / 2
        ns.append(n)
        xs.append('%.2f'%xn)
        ys.append('%.8f'%yn)

四阶龙格库塔

#四阶龙格库塔
def Runge_kutta(x0,y0,rlimX,h,fStr):
    xn = x0
    yn = y0
    n = 0
    ns = [n]
    xs = [xn]
    ys = ['%.8f'%yn]
    while xn < rlimX:
        n += 1
        xn += h
        k1 = globals()[fStr](xn,yn)
        k2 = globals()[fStr](xn + (h / 2),yn + (h * k1) / 2)
        k3 = globals()[fStr](xn + (h / 2),yn + (h * k2) / 2)
        k4 = globals()[fStr](xn + h,yn + h * k3)
        yn = yn + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
        ns.append(n)
        xs.append('%.2f'%xn)
        ys.append('%.8f'%yn)
    return (ns,xs,ys)

演示示例

def ft(x,y):
    return x**2 + x -y
xs = np.linspace(0,5,51)
qt_E1 = EulerExplicit(0,0,0.1,"ft")
qt_E2 = EulerImplicit(0,0,0.1,"ft")
qt_E3 = EulerImproved(0,0,5,0.1,"ft")
qt_R = Runge_kutta(0,0,4.9,0.1,"ft")
plt.figure( figsize=(24,16) )  
print(qt_R[2])
for i in range(51):
    qt_E1[2][i] = float(qt_E1[2][i])
    qt_E2[2][i] = float(qt_E2[2][i])
    qt_E3[2][i] = float(qt_E3[2][i])
    qt_R[2][i] = float(qt_R[2][i])
plt.plot(xs,qt_E1[2],label = "显式欧拉")
plt.plot(xs,qt_E2[2],label = "隐式欧拉")
plt.plot(xs,qt_E3[2],label = "欧拉改进")
plt.plot(xs,qt_R[2],label = "四阶龙格库塔")
plt.legend()
plt.show()

结果如下：