myeuler_vs_ode45.m

比较欧拉公式与matlab的ode方法解微分方程的不同

function yy=myeuler_vs_ode45(fun,xspan,y0,x_step)
% 利用改进欧拉单步迭代法自己编写一阶微分方程的
% 求解函数并与ODE45方法进行比较
% 调用方法：myeuler_vs_ode45(fun,xspan,y0,x_step)
% 参数说明：fun：一阶微分方程描述函数文件名，
%                  以字符串形式给出
%                配合本例老师已经给出了四个微分方程文件，分别
%                是myfun.m, mufun1.m, myfun2.m和myfun3.m
%           xspan:二维行向量，用于指定所感性趣的求解自变量区间
%           y0: 在上述区间左端点处待求解函数的初值
%           x_step: 自编的单步欧拉迭代法的计算迭代步长
% 
scr=get(0,'screensize');     %获得当前显示屏幕的分辨率尺寸
h=clf;                      %清除绘图窗口并获得其图形窗口句柄
set(h,'numbertitle','off','menubar','none','name', ...
    'Euler与ODE45             作者：王希连，仅供大家学习参考', ...
    'position',scr+[0 31 0 -60]);
% 清除标题栏中的图形窗口序号，菜单栏，标题栏中
% 显示'Euler与ODE45'
xx=xspan(1):x_step:xspan(2);
% 根据求解区间和迭代步长计算所有的自变量点
m=length(xx);
yy=zeros(m,1);   % 初始化相应的数值解向量
yy(1)=y0;        
% 25-35行是利用教材上46页（2.5.22式）编写迭代计算数值解向量程序
% 请注意函数feval的使用，并比较它和前面所学过的subs的差别
for n=1:m-1
    x=xx(n);
    y=yy(n);
    k1=feval(fun,x,y);
    x=xx(n+1);
    y=yy(n)+x_step*k1;
    k2=feval(fun,x,y);
    yy(n+1)=yy(n)+x_step/2*(k1+k2);
end
hold on;
[xxx,yyy]=ode45(fun,xspan,y0);    %用ode45函数求解同一个微分方程
% 绘制欧拉迭代法的数值解曲线，我们取出该曲线的图形句柄用于后续控制操作
ymin=min(min(yy),min(yyy));
ymax=max(max(yy),max(yyy));
h1=plot(xx,yy,'Marker','o','Markersize',10,'color',[0.1 0.85 0.1]);
if (abs(max(yyy)-max(yy))<10) &(abs(min(yyy)-min(yy))<10)
    axis([xspan ymin ymax]);      % 固定绘图窗口大小
end
% 给出关于该曲线的文字说明，同样取出它的图形句柄用于后续控制操作
s=strcat('利用步长为',num2str(x_step),'的改进单步欧拉迭代法求解');
h2=text(xspan(1)+diff(xspan)/60,(max(yyy)+min(yyy))/2,s,'FontName', ...
    '华文行楷','Fontsize',26);
pause(2);          %暂停程序执行2秒钟
set(h1,'visible','off');        %隐藏欧拉迭代法的数值解曲线
set(h2,'visible','off');        %隐藏欧拉迭代法的数值解曲线文字说明
% 绘制ODE45的数值解曲线，同样取出该曲线的图形句柄用于后续控制操作
h3=plot(xxx,yyy,'Marker','*','Markersize',10,'color',[0.85 0.2 0.1]);
% 给出关于该曲线的文字说明，同样取出它的图形句柄用于后续控制操作
h4=text(xspan(1)+diff(xspan)/60,(max(yyy)+min(yyy))/2, ...
    '用ODE45方法求解微分方程', ...
    'FontName','华文行楷','Fontsize',26);
mm=0;
% 在X轴下方给出控制程序结束的文字说明
xlabel('点击绘图区域左右两侧以终止程序', ...
    'FontName','华文行楷', ...
    'Fontsize',20,'color',[.2 .1 .95]);
% 使用无限循环语句控制两条曲线的及其相应的文字
% 说明的交替显示，以资比较
while 1
    [gx,gy]=ginput(1);     %获得图形窗口内鼠标的点击位置
    if gx<=xspan(1)|gx>=xspan(2) 
        break;  %如果点击处的x坐标超出求解范围则终止程序
    end
    if(rem(mm,2)==0)
        set(h1,'visible','on');
        set(h2,'visible','on');
        set(h4,'visible','off');
        set(h3,'visible','off');
    else
       set(h3,'visible','on');
       set(h4,'visible','on');
       set(h2,'visible','off');
       set(h1,'visible','off');
   end
    mm=mm+1;
end
close;