3种欧拉法的比较

一、实验内容

常微分方程初值问题的向前欧拉方法、向后欧拉方法、改进欧拉方法数值方法，并以定解问题

定解问题 y'=y-2x/y , y(0)=1

为例，

二、目的与要求

将求解区间进行n=10、20、40、80等分求解

1. 输出剖分n=80时的数值解与精确解图像；
2. 画出误差函数与剖分步长的对数图像；
3. 列表表示n=20时的最大误差、L2误差、L1误差；
4. 写出Matlab代码（包含变量及语句等的注释）;

三、实验步骤（过程）

1. 数值结果

1.1向前欧拉：

 

1.2向后欧拉：

 

 

1.3改进欧拉：

 

 

 

 

代码：

1.主函数

clc;clear;

close all;

%% I.V.P参数初始化

y_exa=@(x) sqrt(2.*x+1);     %原函数

rhs_y=@(x,y) y-2*x/y;       %导函数，方程

y0=1;                       %函数初值

sec=[0,1];                  %解区间

flag=1;                     %{1,2,3}分别对应向前欧拉，向后欧拉，改进欧拉

%% 输出数值解;20剖分下最大误差，L2,L1;四种剖分下的平均误差；误差汇总

[u,E,merr,err]=Euler_3(flag,rhs_y,y_exa,sec,y0);

disp('  最大误差  L2   L1');

E

2.子函数

function [u,E2,perr,err] = Euler_3(flag,rhs_y,y_exa,sec,y0 )

%[数值解，E2,误差均值，误差汇总]={导函数，原函数，解区间，左端初值}

%在四个剖分下计算给定的方程的数值解与和真实解的误差

%by 嫁佛幡

%2022/3/13

%% 初始化参数

t0=sec(1);T=sec(2);                   %解区间

N=[10,20,40,80];            %剖分数

E2=ones(1,3);               %20剖分下最大误差，L2损失函数，L1损失函数

%%

if flag~=1 && flag~=2 && flag~=3

    disp('enter error');

else   

t=cell(4,1);        %四种剖分下的自变量 

u=cell(4,1);        %四种剖分下的数值解

err=cell(4,1);      %四种剖分下的误差汇总

perr=ones(4,1);     %四种剖分下的误差均值

for j=1:4

t{j,1}=linspace(t0,T,N(j)+1);       %自变量赋值

h=(T-t0)/N(j);                      %步长

x=[t0,t0+(1:N(j))*h];               %自变量

y_sol=ones(1,N(j)+1);               %函数值

y_sol(1)=y0;                        %函数初值

if flag==1

mstr='向前欧拉';                    %画图的注释

%y2=y1+h*y1',向前欧拉

for i=1:N(j)

    y_sol(i+1)=y_sol(i)+rhs_y(x(i),y_sol(i))*h;

end

u{j,1}=y_sol;

elseif flag==2

mstr='向后欧拉';                    %画图的注释

%y2=y1+h*y1',向后欧拉

delta=1e-9;                                 %精度

for i=1:N(j)

    y_sol(i+1)=y_sol(i)+rhs_y(x(i),y_sol(i))*h; %向前欧拉

    v0=y_sol(i+1)-1e-8;                         %前点

    v1=y_sol(i+1);                              %后点

    while abs(v1-v0)>delta

        v0=v1;

        v1=y_sol(i)+rhs_y(x(i+1),v0)*h;     %右矩形公式

    end

    y_sol(i+1)=v1;                           %向后欧拉所得值                

end

u{j,1}=y_sol;

elseif flag==3

mstr='改进欧拉';                    %画图的注释

%y2=y1+h*y1',改进欧拉

delta=1e-9;                                 %精度

for i=1:N(j)

    y_sol(i+1)=y_sol(i)+rhs_y(x(i),y_sol(i))*h; %向前欧拉

    v0=y_sol(i+1)-1e-8;                         %前点

    v1=y_sol(i+1);                              %后点

    while abs(v1-v0)>delta

        v0=v1;

        v1=y_sol(i)+(rhs_y(x(i),y_sol(i))+rhs_y(x(i+1),v0))*h/2; %梯形公式 

    end

    y_sol(i+1)=v1;      %改进欧拉所得值               

end

u{j,1}=y_sol;

end

y{j,1}=y_exa(t{j,1});               %解析解

n=length(y{j,1});                    %变量长度                   

L=abs(y{j,1}-u{j,1});               %误差序列

L1=sum(L);                          %误差绝对值和 

L1_mean=mean(L);                    %误差绝对值均值

L2=L*(L')/n;                        %误差绝对值平方均值

e_max=max(L);                       %最大误差绝对值

err{j,1}={L,L1,L1_mean,L2,e_max};   %误差单元赋值

perr(j)=err{j}{1,3};                %误差均值

end

figure(1);

plot(t{4,1},y{4,1},'r-',t{4,1},u{4,1},'b:');        %80剖分下的解析解和数值解图像

title([mstr,'法下80剖分下的解析解和数值解图像']);

xlabel('t');

ylabel('y');

legend('解析解','数值解');

figure(2);

plot(log(1./N),log(perr));                          %误差函数与步长对数的关系

title([mstr,'法下误差函数对数与步长对数的关系']);

xlabel('步长对数');

ylabel('误差均值对数');

E2(1,:)=[err{2}{1,5},err{2}{1,4},err{2}{1,2}];      %20剖分下3个误差赋值

end

end