微分方程数值解法的一些练习
用Euler法解决初值问题
以如下的初值问题为例:
u
′
=
4
t
u
1
/
2
,
0
≤
t
≤
2
u
(
0
)
=
1
\begin{aligned} &u^{'} = 4tu^{1/2},0\leq t \leq 2\\ &u(0)=1 \end{aligned}
u′=4tu1/2,0≤t≤2u(0)=1
该方程的真实解为:
u
(
t
)
=
(
1
+
t
2
)
2
u(t)=(1+t^2)^2
u(t)=(1+t2)2
如下我们进行数值求解,将区间[0,2]进行分划,步长h=0.1则:
f
(
t
n
,
u
n
)
=
4
n
h
u
n
1
/
2
f(t_{n},u_{n})=4nhu_n^{1/2}
f(tn,un)=4nhun1/2
function [y] = f(u,t)
%u的导函数
y = 4*t*sqrt(u);
end
显式欧拉法
公式如下:
u
n
+
1
=
u
n
+
h
f
(
t
n
,
u
n
)
u_{n+1} =u_{n}+hf(t_{n},u_{n})
un+1=un+hf(tn,un)
代码如下:
function [x,y] = Euler(f,a,b,h,y0)
%Euler算法
%f为u的导函数,[a,b]为区间,h为步长,y0为初值
N=(b-a)/h;
x=a:h:b;
y=zeros(1,N+1);
y(1)=y0;
for i=1:N
y(i+1)=y(i)+h*f(y(i),(i-1)*h);
end
end
隐式欧拉法
公式如下:
u
n
+
1
=
u
n
+
h
f
(
t
n
+
1
,
u
n
+
1
)
u_{n+1} =u_{n}+hf(t_{n+1},u_{n+1})
un+1=un+hf(tn+1,un+1)
我们采用预估-校正的方法进行求解,以显格式作为预优算式( P )也即如下形式:
P
:
u
n
+
1
[
0
]
=
u
n
[
M
]
+
h
f
n
[
M
]
E
:
f
n
+
1
[
m
]
=
f
(
t
n
+
1
,
u
n
+
k
[
m
]
)
C
:
u
n
+
1
[
m
+
1
]
=
u
n
M
+
h
f
n
+
1
[
M
]
,
0
≤
m
≤
M
−
1
E
:
f
n
+
1
[
M
]
=
f
(
t
n
+
1
,
u
n
+
k
[
M
]
)
\begin{aligned} &P:u_{n+1}^{[0]} =u_{n}^{[M]}+hf_{n}^{[M]}\\ &E:f_{n+1}^{[m]}=f(t_{n+1},u_{n+k}^{[m]})\\ &C:u_{n+1}^{[m+1]}=u_n^{M}+hf_{n+1}^{[M]},0\leq m \leq M-1\\ &E:f_{n+1}^{[M]}=f(t_{n+1},u_{n+k}^{[M]})\\ \end{aligned}
P:un+1[0]=un[M]+hfn[M]E:fn+1[m]=f(tn+1,un+k[m])C:un+1[m+1]=unM+hfn+1[M],0≤m≤M−1E:fn+1[M]=f(tn+1,un+k[M])
代码如下:
function [x,y] = imp_Euler(f,a,b,h,y0,m)
%impEuler算法
% f为u的导函数,[a,b]为区间,h为步长,y0为初值,m为教证次数
N=(b-a)/h;
x=a:h:b;
y=zeros(1,N+1);
y(1)=y0;
for i=1:N
%P
if i==1
y(i+1)=y(i)+h*f(y(i),(i-1)*h);
else
y(i+1)=y(i)+h*E;
end
for j=1:m
%E
E=f(y(i+1),i*h);
%C
y(i+1)=y(i)+h*E;
end
%E
E=f(y(i+1),i*h);
end
end
改进Euler法(梯形公式)
公式如下:
u
n
+
1
=
u
n
+
h
2
[
f
(
t
n
,
u
n
)
+
f
(
t
n
+
1
,
u
n
+
1
)
]
u_{n+1} =u_{n}+\frac{h}2[f(t_{n},u_{n})+f(t_{n+1},u_{n+1})]
un+1=un+2h[f(tn,un)+f(tn+1,un+1)]
同样利用预估-校正算法求解
代码如下:
function [x,y] = improve(f,a,b,h,y0,m)
%改近Euler算法
%f为u的导函数,[a,b]为区间,h为步长,y0为初值,m为教证次数
N=(b-a)/h;
x=a:h:b;
y=zeros(1,N+1);
y(1)=y0;
for i=1:N
%P
if i==1
y(i+1)=y(i)+0.5*h*f(y(i),(i-1)*h);
else
y(i+1)=y(i)+0.5*h*E;
end
for j=1:m
%E
E=f(y(i+1),i*h);
%C
y(i+1)=y(i)+0.5*h*E+0.5*h*f(y(i),(i-1)*h);
end
%E
E=f(y(i+1),i*h);
end
end
main.m如下:
[x,y1]=Euler(@f,0,2,0.1,1);
[~,y2]=imp_Euler(@f,0,2,0.1,1,2);
[~,y3]=improve(@f,0,2,0.1,1,2);
g=@(z)(1+z.^2).^2;
z=x;
y=g(z);
figure(1)
plot(x,y1);
hold on
plot(x,y2)
plot(x,y3)
plot(x,y)
scatter(x,y1)
scatter(x,y2)
scatter(x,y3)
scatter(x,y,'*')
legend('E','impE','improve','gt','Ep','impEp','imporvep','gtp')
title('Euler法相关(隐式采用P(EC)^mE)')
hold off
结果图
参考文献
[1]李荣华,刘播. 微分方程数值解法(第四版).北京:高等教育出版社,2009.