Euler法解微分方程

用Euler法解决初值问题

以如下的初值问题为例：
$$\begin{array}{rl}& u^{{}^{\prime }}=4t{u}^{1/2},0\le t\le 2\\ & u(0)=1\end{array}$$
该方程的真实解为：
$$u(t)=(1+t^2{)}^2$$
如下我们进行数值求解,将区间[0,2]进行分划，步长h=0.1则：
$$f(t_n,u_n)=4nh{u}_n^{1/2}$$

function [y] = f(u,t)
%u的导函数
y = 4*t*sqrt(u);
end

显式欧拉法

公式如下：
$$u_{n+1}=u_n+hf(t_n,u_n)$$
代码如下：

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+hf(t_{n+1},u_{n+1})$$
我们采用预估-校正的方法进行求解，以显格式作为预优算式( P )也即如下形式：
$$\begin{array}{rl}& 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\le m\le M-1\\ & E:f_{n+1}^{[M]}=f(t_{n+1},u_{n+k}^{[M]})\end{array}$$
代码如下：

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+\frac{h}{2}[f(t_n,u_n)+f(t_{n+1},u_{n+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.