# 微分方程数值解法

## 改进Euler法

### 局部误差分析

\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}

### 整体误差

$$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)$$

$$|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)}$$

### 稳定性

$$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+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)]$$

\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}

### 代码

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]))