欧拉数值方法(lesson 2)

一种一阶数值方法，用以对给定初值的常微分方程求解

对于一阶常系数微分方程 $y'=f(x,y)$

理论推导:

二元函数 $f(x,y)$可以看做y的斜率
给定初始值 $y(a) = b$ ,即 $(a,b)$, 步长 $h$,可以求出下一临近点的位置
用公式表示:

$x_{n+1}=x_{n}+h$

$y_{n+1}=y_{n}+h*An$

$An=f(x_{n},y_{n})$

不断进行迭代就可以汇出 $y$的图像
此时的图像误差较大
当 $y''<0$时欧拉数值法的结果会偏大
当 $y''>0$时欧拉数值法的结果会偏小

优化欧拉数值法

1. 减小步长 $h$
2. 采用更高阶(这里以四阶为例)

$x_{n+1}=x_{n}+h$

$y_{n+1}=y_{n}+h*(An+2Bn+2Cn+Dn)*1/6$

$An=f(x_{n},y_{n})$

$Bn=f(x_{n}+h,h*An)$

$Cn=f(x_{n}+2h,h*(An+Bn))$

$Dn=f(x_{n}+3h,h*(An+Bn+Cn))$

import matplotlib.pyplot as plt

#y' = x^2 - y^2
#f(x, y)=x^2 - y^2
def An(x, y):
    return x**2 - y**2
def Bn(x, y, h):
    return (x + h)**2 - (y + h*An(x, y))**2
def Cn(x, y, h):
    return (x + 2*h)**2 - (y + h*(An(x, y) + Bn(x, y, h)))**2
def Dn(x, y, h):
    return (x + 3*h)**2 - (y + h*(An(x, y) + Bn(x, y, h) + Cn(x, y, h)))**2
def One():
    X = []
    Y = []
    x0, y0 = 0, 1
    iteration = 10000
    h = 0.0002
    for i in range(iteration):
        x1 = x0 + h
        y1 = y0 + h*An(x0, y0)
        X.append(x1)
        Y.append(y1)
        x0, y0 = x1, y1
    plt.plot(X, Y)
    plt.title('One')
    plt.show()
def Two():
    X = []
    Y = []
    x0, y0 = 0, 1
    iteration = 10000
    h = 0.0002
    for i in range(iteration):
        x1 = x0 + h
        y1 = y0 + h*(An(x0, y0) + Bn(x0, y0, h))*.5
        X.append(x1)
        Y.append(y1)
        x0, y0 = x1, y1
    plt.plot(X, Y)
    plt.title('Two')
    plt.show()

def Four():
    X, Y = [], []
    x0, y0 = 0, 1
    iteration = 10000
    h = .0002
    for i in range(iteration):
        x1 = x0 + h
        y1 = y0 + h*(An(x0,y0) + 2*Bn(x0,y0,h) 
                    + 2*Cn(x0,y0,h) + Dn(x0,y0,h))*(1/6)
        X.append(x1)
        Y.append(y1)
        x0, y0 = x1, y1
    plt.plot(X, Y)
    plt.title('Four')
    plt.show()

if __name__ == '__main__' :
    One();
    Two();
    Four();