欧拉近似方法求常微分方程.doc

《欧拉近似方法求常微分方程.doc》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、欧拉近似方法求常微分方程朱翼1、编程实现以下科学计算算法，并举一例应用之。“欧拉近似方法求常微分方程”算法说明：欧拉法是简单有效的常微分方程数值解法，欧拉法有多种形式的算法，其中简单欧拉法是一种单步递推算法。其基本原理为对简单的一阶方程的初值问题：y’=f(x,y)其中y(x0)=y0欧拉法等同于将函数微分转换为数值微分，由欧拉公式可得yn+1=yn+hf(xn,yn)程序代码：function[tout,yout]=myeuler(ypfun,t0,tfinal,y0,tol,trace)%初始化pow=1/3;ifnargin<5,tol=1.e-

2、3;endifnargin<6,trace=0;endt=t0;hmax=(tfinal-t)/16;h=hmax/8;y=y0(:);chunk=128;tout=zeros(chunk,1);yout=zeros(chunk,length(y));k=1;tout(k)=t;yout(k,:)=y.';iftrace%绘图clc,t,h,yendwhile(tt)%主循环ift+h>tfinal,h=tfinal-t;end%Computetheslopesf=feval(ypfun,t,y);f=f(:);%估计误差并

3、设定可接受误差delta=norm(h*f,'inf');tau=tol*max(norm(y,'inf'),1.0);%当误差可接受时重写解ifdelta<=taut=t+h;y=y+h*f;k=k+1;ifk>length(tout)tout=[tout;zeros(chunk,1)];yout=[yout;zeros(chunk,length(y))];endtout(k)=t;yout(k,:)=y.';endiftracehome,t,h,yend%Updatethestepsizeifdelta~=0.0h=min(hmax,0.9*h*(

4、tau/delta)^pow);endendif(t

5、k,length(y));k=1;tout(k)=t;yout(k,:)=y.';NYtrace=1N输出t,h,yNttYYt+h>tfinalh=tfinal-tNf=feval(ypfun,t,y);f=f(:);delta=norm(h*f,'inf');tau=tol*max(norm(y,'inf')1.0);tau=tol*max(norm(y,'inf')1.0);Ydelta<=tauNt=t+h;y=y+h*f;k=k+l;k>length(tout)tout=[tout;zeros(chunk,1)];yo

6、ut=[yout;zeros(chunk,length(y))];tout(k)=t;yout(k,:)=y,';Ytrace=1输出home,t,h,yNYdelta~=0.0h=min(hmax,0.9*h*(tau/delta)^pow))NYt

7、=x+exp(-x)。最后在同一图窗中分别画出两图。程序代码：f．mfunctionf=f(x,y)f=-y+x+1;>>[x,y]=myeuler('f',0,1,1);%利用程序求解方程>>y1=x+exp(-x);%方程f=-y+x+1的精确解>>plot(x,y,'-b',x,y1,'-r')%在同一图窗将欧拉法解和精确解的图画出>>legend('欧拉法','精确解')例题流程图：f=f(x,y)f=-y+x+1;开始调用函数myeuler，[x,y]=myeuler(‘f’,0,1,1);求出结果y1=x+exp(-x)调用函数plot绘图

8、，比较结束2、编程解决以下科学计算问题（1）解题思路：◆建模：材料力学中从弯矩求转角要经过一次