Featured image of post 欧拉方法、隐式欧拉方法

欧拉方法、隐式欧拉方法

求解如下常微分方程:

$$ \begin{aligned} \begin{cases} \frac{du}{dt}=-\frac{u}{t},1\leq t\leq 2, \\\ u(1)=1 \end{cases} \end{aligned} $$

精确解

将原方程化为$tdu+udt=0$,则有$d(ut)=0$,解得$ut=C$($C$为常数),代入初始条件得$C=1$,从而该方程的精确解为:

$$ \begin{aligned} u=\frac{1}{t},(1\leq t \leq2). \end{aligned} $$

欧拉方法

代入欧拉格式得:

$$ \begin{aligned} u_{i+1}=u_{i}+hf(t_i,u_i)=u_i+h(-\frac{u_i}{t_i}) \end{aligned} $$

隐式欧拉方法

由隐式欧拉格式得:

$$ \begin{aligned} u_{i+1}=u_{i}+hf(t_{i+1},u_{i+1})=u_i+h(-\frac{u_{i+1}}{t_{i+1}}), \end{aligned} $$

移项化简可得:

$$ \begin{aligned} u_{i+1}=\frac{t_{i+1}u_i}{t_{i+1}+h} \end{aligned} $$

程序

根据上述推导,用python编写程序,代码如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
# implict euler method
import numpy as np
import matplotlib.pyplot as plt

# the right term of the ODE
def f(t, u):
    f = -u/t
    return f

# the exact solution of ODE 
def fexact(t):
    fexact = 1/t
    return fexact

N = 100
t_n = 2.0
dt = (t_n - 1.0) / N
t = np.arange(1.0, t_n + dt, dt)
u_euler = np.arange(1.0, t_n + dt, dt)
u = np.arange(1.0, t_n + dt, dt)
u_true = np.arange(1.0, t_n + dt, dt)

i = 0
while i < N:
    t[i+1] = t[i] + dt
    u_euler[i+1] = u_euler[i] + dt * f(t[i], u_euler[i])
    u[i+1] = (u[i] * t[i+1])/(t[i+1] + dt)
    u_true[i+1] = fexact(t[i+1])
    i = i + 1

err_euler = max(abs(u_euler - u_true))
err_implict_euler = max(abs(u - u_true))
print("The error of euler method: ",err_euler)
print("The error of implict euler method: ",err_implict_euler)

# begin drawing
plt.title('Result')
plt.plot(t, u_euler, color='green', label='euler')
plt.plot(t, u, color='blue', label='implict euler')
plt.plot(t, u_true, color='red', label='exact')
plt.legend() # show the legend

plt.xlabel('t')
plt.ylabel('u')
plt.show()

结果分析

当取$h=0.01$时,此时欧拉方法的误差为0.02631578947368396,隐式欧拉方法的误差为0.023809523809523836,结果如下图所示:

result.png

当取不同$h$,得到的误差如下表所示:

$h$欧拉方法隐式欧拉方法
$\frac{1}{2}$0.166666666666666630.09999999999999998
$\frac{1}{4}$0.07142857142857140.05555555555555558
$\frac{1}{8}$0.0333333333333332150.02941176470588236
$\frac{1}{16}$0.016129032258064670.015151515151515138
$\frac{1}{32}$0.007936507936508130.007692307692307665
$\frac{1}{64}$0.00393700787401551930.003875968992248069