欧拉数值方法(lesson 2)

欧拉数值方法(lesson 2)

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

对于一阶常系数微分方程$y^{\prime }=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\ast An$

$An=f(x_n,y_n)$

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

优化欧拉数值法

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

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

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

$An=f(x_n,y_n)$

$Bn=f(x_n+h,h\ast An)$

$Cn=f(x_n+2h,h\ast (An+Bn))$

$Dn=f(x_n+3h,h\ast (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();