matlab 常微分方程数值解法PPT课件

.,科学计算与MATLAB,主讲：唐建国中南大学材料科学与工程学院2012,.,第七讲常微分方程数值解法,.,内容提要,引言欧拉近似方法龙格-库塔(R-K)方法MATLAB的常微分方程函数小结,.,1、引言,物理学所研究的各种物质运动中，有许多物质运动的过程是用常微分方程来描述的。,例如,质点的加速运动，简谐振动等。,简单问题可以求得解析解，多数实际问题靠数值求解。,.,一阶常微分方程(ODE)初值问题：,数值解法就是求y(x)在某些分立的节点xn上的近似值yn，用以近似y(xn),ODE：OrdinaryDifferentialEquation,.,2.1简单欧拉(L.Euler,1707-1783)方法。,欧拉数值算法就是由初值通过递推求解，递推求解就是从初值开始，后一个函数值由前一个函数值得到。关键是构造递推公式。,2、欧拉近似方法,.,欧拉数值算法递推公式构造,2.1.1差分法,差分法就是用差商近似代替微商，即,代入微分方程得到：,.,对于等间隔节点,可以得到：,.,在xn节点上，微分方程可以写为,作如下近似：,则得到欧拉解法递推公式的一般形式：,.,具体求解过程为：,.,简单欧拉方法程序functionoutx,outy=MyEuler(fun,x0,xt,y0,PointNum)%MyEuler用前向差分的欧拉方法解微分方程%fun表示f（x，y）%x0,xt表示自变量的初值和终值%y0表示函数在x0处的值,其可以为向量形式%PointNum表示自变量在x0,xt上取的点数ifnargin5|PointNumdsolve(Dy+3*x*y=x*exp(-x2)ans=(1/3*exp(-x*(x-3*t)+C1)*exp(-3*x*t)dsolve(Dy+3*x*y=x*exp(-x2),x)ans=exp(-x2)+exp(-3/2*x2)*C1,.,例题：用MATLAB的符号解法,求解常微分方程的特解：,dsolve(x*Dy+2*y-exp(x)=0,y(1)=2*exp(1),x)ans=(exp(x)*x-exp(x)+2*exp(1)/x2,.,例题：采用ODE45求解描述某刚性问题的微分方程：,functiondy=rigid(t,y)dy=zeros(3,1);%行向量dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);,options=odeset(RelTol,1e-4,AbsTol,1e-41e-41e-5);T,Y=ode45(rigid,012,011,options);plot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.)legend(y1,y2,y3),.,.,例题：用MATLAB函数ode23,ode45,ode113,求解多阶常微分方程：,.,functiondy=myfun03(x,y)dy=zeros(3,1)%初始化变量dydy(1)=y(2);%dy(1)表示y的一阶导数,其等于y的第二列值dy(2)=y(3);%dy(2)表示y的二阶导数dy(3)=2*y(3)/x3+3*y(2)/x3+3*exp(2*x)/x3%dy(3)表示y的三阶导数,%用ode23ode45ode113解多阶微分方程clear,clcx23,y23=ode23(myfun03,1,10,11030);x45,y45=ode45(myfun03,1,10,11030);x113,y113=ode113(myfun03,1,10,11030);figure(1)%第一幅图plot(x23,y23(:,1),*r,x45,y45(:,1),ob,x113,y113(:,1),+g)%作出各种函数所得结果legend(ode23解,ode45解,ode113解)title(ODE函数求解结果)figure(2)plot(x45,y45)%以ode45为例作出函数以及其各阶导数图legend(y,y一阶导数,y两阶导数)title(y,y一阶导数,y二阶导数函数图),.,.,例题：MATLAB在导热计算,传热过程涉及面广，数学模型复杂。计算过程中涉及到许多运算方法。导热虽是容易做数学处理的一种热量传递方式，但其过程往往涉及常微分、偏微分方程、线性（非线性）方程组的求解。,.,有一外径为4cm，内径为1.5cm，载有电流密度I为5000A/cm2的内冷钢质导体。导体单位时间发出的热量等于流体同时带走的热量，导体内壁面的温度为70。假定外壁面完全绝热。试确定，导体内部的温度分布；已知钢的导热系数k=0.38Kw/mK，电导率=21011km。,解：这是圆柱坐标中常物性一维稳态导热问题，结合本题具体条件导热微分方程式可简化为：,.,结合边界条件，可得这一导热问题的数学描述为：,此常微分方程的分析解，可调用MATLAB符号工具箱中的dsolve函数来实现。在命令窗口执行下面的代码：cleard_equat=D2t+Dt/r+131579=0;condition=（t0.0075）=70,D（t0.02）=0;%边界条件t=dsolve（d_equat,condition,r）程序执行结果：,.,程序执行结果：,即求出温度分布方程为:,.,工程上遇到的导热问题，往往由于物体的几何形状复杂或边界条件难以描述，无法求出分析解，此时可采用数值方法进行求解。常微分方程（ODE）包括初值问题和边值问题两种，初值问题ODE的数值解法常调用ode45（）和ode23（）函数实现。在MATLAB编辑器中编写函数BZ.m，存盘。functionBZclearallclca=0.0075;b=0.02;%r值的范围solinit=bvpini（tlinspace（a,b,10）,7080）；%用初始值对解初始化sol=bvp4c（ODEfun,BCfun,solinit）；%解ODE方程r=0.0075:0.00125:0.02t=deva（lsol,r）functiondtdr=ODEfun（r,t）%定义ODE方程函数dtdr=（t2）;-131580-1/r