微分方程数值解法

改进Euler法

局部误差分析

将$t\in[t_n,t_{n+1}]$表示成$t=t_n+\tau h,0\leq\tau\leq 1$. 由线性插值的余项公式，我们有
$$
\begin{aligned}f(t,u(t))&=u'(t)=u'(t_n+\tau h)\\&=u'(t_n)+u''(t_n)(\tau h)+\frac{h^2}{2}u'''(t+\theta h)\omega_2(t)\\&=u'(t_n)+\tau[u'(t_{n+1})-u'(t_{n})]+\frac{h^2}{2}\tau(\tau-1)u'''(t_n+\theta h),(0\leq\theta\leq1)\end{aligned}
$$

插值余项及其估计

于是
$$
\begin{aligned}\int_{t_n}^{t_{n+1}}f(t,u(t))dt&=\int_0^1f(2h,u(2h))hd\tau\\&=\int_0^1\left(u'(t_n)+\tau[u'(t_{n+1})-u'(t_{n})]\right)d\tau+\frac{h^3}{2}\int_0^1\tau(\tau-1)u'''(t_n+\theta h)d\tau\\&=\frac{h}{2}[u'(t_n)+u'(t_{n+1})]-\frac{h^3}{12}u'''(\zeta), \zeta\in(t_n,t_{n+1}).\end{aligned}
$$
改进Euler法$u(t_1)=u_0+\int_{t_0}^{t_1}f(\tau,u(\tau))d\tau$选择了梯形公式近似$u_1=u_0+\frac{h}{2}[f(t_0,u_0)+f(t_1,u_1)]$. 所以其局部截断误差为$R_n^{(1)}=-\frac{h^3}{12}u'''(\zeta)$, 其阶为$O(h^3)$.

整体误差

$$
u(t_n+1)=u(t_n)+\int_{t_0}^{t_1}f(\tau,u(\tau))d\tau\;(t_0=0).
$$

$$
u_{n+1}=u_n+\frac{h}{2}\left[f(t_n,u_n)+f(t_{n+1},u_{n+1})\right]
$$

两个式子做差
$$
|e_{n+1}|=|e_n|-\frac{h^3}{12}u'''(\zeta)
$$
记$R^{(1)}=\max R^{(1)}_n$, 其中$t_0+nh<T$
$$
|e_{n+1}|\leq|e_n|-R^{(1)}\leq e_0+(n+1)R^{(1)}=e_0+\frac{T-t_0}{h}R^{(1)}+R^{(1)}
$$
所以改进Euler格式的整体误差是$O(h^2)$

稳定性

$$
u_{n+1}=u_{n}+\frac{h}{2}[f(t_n,u_n)+f(t_{n+1},u_{n+1})]
$$

$$
v_{n+1}=v_n+\frac{h}{2}[f(t_n,v_n)+f(t_{n+1},v_{n+1})]
$$

两式相减并令$e_n=u_n-v_n$, 得
$$
e_{n+1}=e_{n}+\frac{h}{2}[\left(f(t_n,u_n)-f(t_n,v_n)\right)+\left(f(t_{n+1},v_{n+1})-f(t_{n+1},v_{n+1})\right)]
$$
因满足Lipschitz条件从而
$$
\begin{aligned}|e_{n+1}|&\leq |e_n|+\frac{h}{2}[L|e_n|+L|e_{n+1}|]\\&\leq\left(\frac{1+Lh}{1-Lh}\right)^n|e_0|\leq e^{LT}|e_0|\end{aligned}
$$
因$nh\leq T$这说明$e_n$连续依赖初始误差$e_0$. 即Euler方法是稳定的.

代码

def compute_fx(x,y):
    return y-2*x/y

def Modified_Euler(x_initial,y_initial,h1,N1):
    n=0
    y1=[]
    while n<N1:
        n+=1
        x1=x_initial+h1
        y1_pre=y_initial+h1*compute_fx(x_initial,y_initial)
        y1_modify=y_initial+(compute_fx(x_initial, y_initial)+compute_fx(x1, y1_pre))*h1/2
        y1.append(y1_modify)
        x_initial=x1
        y_initial=y1_modify
    return y1

x0=y0=h=N=0
x0,y0,h,N=map(float,input().strip().split())
result=[]
result=Modified_Euler(x0, y0, h, N)
for i in range(1,int(N)+1):
    X=x0+h*i
    print("{x:.3f}  -->  {y:.5f}".format(x=X, y=result[i-1]))