科学工程计算与matlab编程5-2讲解.ppt

5.2 微分方程问题的数值解法; 微分方程求解的误差与步长问题：;舍入误差;提高数值解精度的方法之一： 减小步长h的值；; function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) % fun 表示f(x,y); x0,xt:自变量的初值和终值; %y0:函数在x0处的值,其可以为向量形式; %PointNum表示自变量在[x0,xt]上取的点数 if nargin==4 | PointNum=0 %PointNum 默认值为100 PointNum=100; end if nargin==3 % y0默认值为0 y0=0; PointNum=100; end;h=(xt-x0)/PointNum; %计算步长h x=x0+[0:PointNum]*h; %自变量数组 y= zeros(size(x));％与自变量点相应的函数值 y(1)=y0; %初值 for k = 1:PointNum y(k+1)= y(k)+h*feval(fun,x(k),y(k)); %对于所取的点x迭代计算y值 end outy=y; outx=x; plot(x,y) %画出方程解的函数图;例：求;y=dsolve(Dy=y+sin(t),y(0)=1); %该常微分方程的解析解 for k=1:33 t(k)=x11(k); y2(k)=subs(y,t(k)); %求其对应点的离散解 end plot(x1,y1,+b,x11,y11,og,x11,y2,*r) legend(h=0.4的欧拉法解,h=0.2的欧拉解,符号解);;改进的Euler算法;function [Xout,Yout]=MyEulerPro(fun,x0,xt,y0,PointNumber) %MyEulerPro 用改进的欧拉法解微分方程 if nargin==4 | PointNumber=0 %PointNumer默认值为100 PointNumer=100; end if nargin==3 %y0默认值为0 y0=0; PointNumer=100; end;h=(xt-x0)/PointNumber; %计算所取的两离 散点之间的距离 x=x0+[0:PointNumber]*h; %表示出离散的自变量x,初始点为x0 y=zeros(size(x)); ％输入与x相应的函数值向量 y(1)=y0;％初始值 for k=1:PointNumber %迭代计算过程 z1=y(k)+h*feval(fun,x(k),y(k)); z2=y(k)+h*feval(fun,x(k+1),z1); y(k+1)=1/2*(z1+z2); end Xout=x;Yout=y;;例：对方程 应用改进欧拉法、简单欧拉方法 以及微分方程符号求解方法分别求解并比较结果； f1=@(x,y)sin(x)+y; [x3,y3]=MyEulerPro(f1,0,2*pi,1,128); [x,y1]=MyEuler(f1,0,2*pi,1,128);%欧拉法所的的解 y=dsolve(Dy=y+sin(t),y(0)=1);%该常微分方程的符号解 for k=1:129 t(k)=x(k); y2(k)=subs(y,t(k));%求解析解对应点的离散解 end plot(x,y1,-b,x3,y3,og,x,y2,*r) legend(简单欧拉法解,改进欧拉解,符号解); ;5.2.2 四阶定步长Runge-Kutta算法及Matlab实现;例如：p＝4，即;Matlab 实现： function [tout,yout]=rk_4(odefile,tspan,y0) ％y0初值列向量 t0=tspan(1); th=tspan(2); if length(tspan)=3, h=tspan(3); % tspan=[t0,th,h] else, h=tspan(2)-tspan(1); th=tspan(end); end ％等间距数组 tout=[t0:h:th]; yout=[]; for t=tout k1=h*eval对应的是数组([odefil