后退 Euler 求解非线性 ODE

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

后退 Euler 方法

$$u_{n+1}=u_{n}+hf(t_{n+1},u_{n+1}).$$

线性 ODE 例子（做对比）

$$\left\{\begin{aligned} &\frac{du}{dt}=t^2+t-u ,\quad 0 < t \leqslant 1, \\ &u(0)=0. \end{aligned} \right.$$

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

$$u_{n+1}=u_n+h(t_{n+1}^2+t_{n+1}-u_{n+1}).$$ $$u_{n+1}=\frac{1}{1+h}(u_n+ht_{n+1}^2+ht_{n+1}).$$

显格式比较容易编程。

非线性 ODE

$$\left\{\begin{aligned} &\frac{\mathrm{d}u}{\mathrm{d}t}=u-\frac{2t}{u} ,\quad 0 < t \leqslant 1, \\ &u(0)=1. \end{aligned} \right.$$

方程的真解：.

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

$$u_{n+1}=u_n+h(u_{n+1}-\frac{2t_{n+1}}{u_{n+1}}),$$

可化为

$$u_{n+1}-u_n-hu_{n+1}+h\frac{2t_{n+1}}{u_{n+1}}=0.$$

每一步要解一个关于 的非线性方程。

牛顿迭代

假设 有近似根 ( )，求 的方法

$$x_{k+1}=x_{k}-\frac{f(x_{k})}{f^{\prime}(x_{k})}, \quad k=0,1, \cdots.$$

对于以上非线性微分方程的离散，记 (为了不混淆)

$$F(X)=X-u_n-hX+\frac{2ht_{n+1}}{X},$$ $$F^{\prime}(X)=1-h-\frac{2ht_{n+1}}{X^2}.$$

数值格式

$$X_{k+1}=X_{k}-\frac{F(X_{k})}{F^{\prime}(X_{k})}, \quad k=0,1,\cdots,N-1，$$ $$X_{k+1}=X_{k}-\frac{X_k-u_n-hX_k+\frac{2ht_{n+1}}{X_k}}{1-h-\frac{2ht_{n+1}}{X_k^2}}, \quad k=0,1,\cdots,N-1.$$

设置终止迭代条件 ( 为允许误差)，最后算下误差。如果有真解的话，算下误差阶。

计算数值解

% 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_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

输出结果

Error =
    0.0702    0.0322    0.0061    0.0030    0.0006

order =
    1.1238    1.0377    1.0103    1.0035

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

$$\left\{\begin{aligned} &\frac{\mathrm{d}u}{\mathrm{d}t}=-100u^3+\sin(u) ,\quad 0 < t \leqslant 1, \\ &u(0)=1. \end{aligned} \right.$$

方程没有真解。

计算数值解

% 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

输出结果