Andy123t
0%

后退 Euler 求解非线性 ODE

发表于 更新于 分类于

本文用后退 Euler 方法求解非线性 ODE，需要用到牛顿迭代，然后用 MATLAB 编程，最后展示了输出结果。

后退 Euler 方法

线性 ODE 例子（做对比）

用后退 Euler (隐式 Euler) 格式离散

显格式比较容易编程。

非线性 ODE

方程的真解：$u(t)=\sqrt{2t+1}$.

用后退Euler (隐式Euler) 格式离散

可化为

每一步要解一个关于 $u_{n+1}$ 的非线性方程。

牛顿迭代

假设 $f(x)=0$ 有近似根 $x_k$ ($f^{\prime}\left(x_{k}\right) \neq 0$ )，求 $x_{k+1}$ 的方法

对于以上非线性微分方程的离散，记 $u_{n+1}=X$ (为了不混淆)

数值格式

设置终止迭代条件 $|X_{k+1}-X_k|<\delta$ ($\delta$ 为允许误差)，最后算下误差。如果有真解的话，算下误差阶。

计算数值解

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
% BackEuler.m
% Backward Euler Method for Nonlinear ODE:
% u'(x)=u-2x/u in [0,1]
% Initial condition: u(0)=1;
% Exact solution: u=sqrt(2*x+1)
clear all; close all;
h=0.01;
x=0:h:1;
N=length(x)-1;
u(1)=1;
NI(1,N)=0;  % Record the number of iterations
for n=1:N
    % Newton iteration
    Xn=u(n);
    Xp=Xn;
    Xprev=0;
    while abs(Xp-Xprev) > eps*abs(Xp)
        Xprev=Xp;
        Xp=Xp-(Xp-h*Xp+2*h*x(n+1)/Xp-u(n))./(1-h-2*h*x(n+1)/(Xp^2));
        NI(n)=NI(n)+1;
    end
    u(n+1)=Xp;
end
ue=sqrt(2*x+1);    % exact solution
plot(x,ue,'b-',x,u,'r+','LineWidth',1)
legend('Exact ','Numerical','location','northwest')
% title('Backward Euler method','fontsize',12)
set(gca,'fontsize',12)
xlabel('t','fontsize',16), ylabel('u','fontsize',16)

% computing error
error=max(abs(u-ue))

% print -dpng -r600 BackEuler.png
% print -depsc2 BackEuler.eps

输出结果

BackEuler.png

计算误差阶

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
46
47
48
49
50
51
52
53
54
% BackEuler_error.m
% Backward Euler Method for Nonlinear ODE:
% u'(x)=u-2x/u in [0,1]
% Initial condition: u(0)=1;
% Exact solution: u=sqrt(2*x+1)
clear all; close all;
Nvec=[10 20 100 200 1000];
%hv=[0.1 0.05 0.01 0.005 0.001];
Error=[];
for k=1:length(Nvec)
N=Nvec(k);
h=1/N;
x=0:h:1;
N=length(x)-1;
u(1)=1;
NI(1,N)=0;  % Record the number of iterations
for n=1:N
    % Newton iteration
    Xn=u(n);
    Xp=Xn;
    Xprev=0;
    while abs(Xp-Xprev) > eps*abs(Xp)
        Xprev=Xp;
        Xp=Xp-(Xp-h*Xp+2*h*x(n+1)/Xp-u(n))./(1-h-2*h*x(n+1)/(Xp^2));
        NI(n)=NI(n)+1;
    end
    u(n+1)=Xp;
end
ue=sqrt(2*x+1);
error=max(abs(u-ue));
Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(log10(Nvec), log10(Nvec.^(-1)), '--')
grid on
%title('Rate of Convergence','fontsize',12)
set(gca,'fontsize',14)
xlabel('log_{10}N','fontsize',14),ylabel('log_{10}Error','fontsize',14)

% add annotation of slope
ax = [0.57 0.53];
ay = [0.68 0.63];
annotation('textarrow',ax,ay,'String','slope = -1 ','fontsize',14)

% computating convergence order
for n=1:length(Nvec)-1     % computating convergence order
    order(n)=-log(Error(n)/Error(n+1))/(log(Nvec(n)/Nvec(n+1)));
end
Error
order

% print -dpng -r600 BackEuler_error.png
% print -depsc2 BackEuler_error.eps

输出结果

1
2
3
4
5
Error =
    0.0702    0.0322    0.0061    0.0030    0.0006

order =
    1.1238    1.0377    1.0103    1.0035

BackEuler_error.png

另一个非线性 ODE 例子 也是可以算的。

方程没有真解。

计算数值解

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
% BackEuler2.m
% Backward Euler Method for Nonlinear ODE:
% u'(t)=-100*u^3+sin(u) in [0,1]
% Initial condition: u(0)=1.
clear all; close all;
h=0.001;
x=0:h:1;
N=length(x)-1;
u(1)=1;
NI(N)=0;  % Record the number of iterations
for n=1:N
    % Newton iteration
    Xn=u(n);
    Xp=Xn;
    Xprev=0;
    while abs(Xp-Xprev) > 1.0e-8*abs(Xp)
        Xprev=Xp;
        Xp=Xp-(Xp+100*h*Xp^3-h*sin(Xp)-u(n))./(1+300*h*Xp^2-h*cos(Xp));
        NI(n)=NI(n)+1;
    end
    u(n+1)=Xp;
end
plot(x,u,'r-.', 'LineWidth',1)
set(gca,'fontsize',12)
xlabel('t','fontsize',16), ylabel('u','fontsize',16)

% print -dpng -r600 BackEuler2.png
% print -depsc2 BackEuler2.eps

输出结果

BackEuler2.png