常微分方程初值问题的数值方法 Python实现

原理：

（1）利用Euler方法和改进的Euler方法求解初值问题
（2）利用Runge-Kutta方法求解初值问题

步骤：

import math

def f1(x, y):
    if x == 0:
        return 0
    else:
        return ((4 * x) / y) - (x * y)
def f2(x, y):
    return x ** 2 - y ** 2


def Euler(a, b, alpha, N):
    h = (b - a) / N
    x = a
    y = alpha
    out = [(x, y, y, y-y)]
    print ("%-10s %-10s %-10s %-10s"%("Xn","Yn","Y(Xn)","Y(Xn)-Yn"))
    print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,y,y-y))
    for n in range(1, N + 1):
        y = y + h * f1(x,y)
        x = x + h
        Y = (4 + 5 * (math.e ** -(x ** 2))) ** 0.5
        out.append((x, y, Y, y-Y))
        print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,Y,y-Y))
    return out

def EulerPlus(a, b, alpha, N):
    h = (b - a) / N
    x = a
    y = alpha
    print ("%-10s %-10s %-10s %-10s"%("Xn","Yn","Y(Xn)","Y(Xn)-Yn"))
    print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,y,y-y))
    out = [(x, y, y, y-y)]
    for n in range(1, N + 1):
        K1 = f1(x, y)
        K2 = f1((x + h), (y + h * K1))
        y = y + h * (K1 + K2) / 2
        x = x + h
        Y = (4 + 5 * (math.e ** -(x ** 2))) ** 0.5
        out.append((x, y, Y, y-Y))
        print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,Y,y-Y))
    return out


def RungeKutta(a, b, alpha, N):
    h = (b - a) / N
    x = a
    y = alpha
    out = [(x,y)]
    print ("%-10s %-10s"%("Xn","Yn"))
    print ("%-10.6f %-10.6f"%(x,y))
    for n in range(1, N + 1):
        K1 = f2(x, y)
        K2 = f2((x + h / 2),(y + (h * K1) / 2))
        K3 = f2((x + h / 2),(y + (h * K2) / 2))
        K4 = f2((x + h),(y + h * K3))
        y = y + h * (K1 + 2 * K2 + 2 * K3 + K4) / 6
        x = x + h
        out.append((x,y))
        print ("%-10.6f %-10.6f"%(x,y))
    return out



Euler(0,2,3,20)
Euler(0,2,3,10)
Euler(0,2,3,5)
EulerPlus(0,2,3,20)
EulerPlus(0,2,3,10)
EulerPlus(0,2,3,5)
RungeKutta(-1,0,0,10)
RungeKutta(-1,0,0,5)