matlab 常微分方程数值解法.ppt

1、A,1,科学计算与MATLAB,主讲：唐建国 中南大学材料科学与工程学院 2012,A,2,第七讲 常微分方程数值解法,A,3,内容提要,引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结,A,4,1、引言,物理学所研究的各种物质运动中，有许多物质运动的过程是用常微分方程来描述的。,例如,质点的加速运动，简谐振动等。,简单问题可以求得解析解，多数实际问题靠数值求解。,A,5,一阶常微分方程(ODE )初值问题 ：,数值解法就是求y(x)在某些分立的节点 xn 上的近似值yn，用以近似y(xn),ODE ：Ordinary Differential Equation

2、,A,6,2.1 简单欧拉(L.Euler, 1707-1783)方法。,欧拉数值算法就是由初值通过递推求解，递推求解就是从初值开始，后一个函数值由前一个函数值得到。关键是构造递推公式。,2、欧拉近似方法,A,7,欧拉数值算法递推公式构造,2.1.1 差分法,差分法就是用差商近似代替微商，即,代入微分方程得到：,A,8,对于等间隔节点,可以得到：,A,9,在xn节点上，微分方程可以写为,作如下近似：,则得到欧拉解法递推公式的一般形式：,A,10,具体求解过程为：,A,11,简单欧拉方法程序 function outx,outy=MyEuler(fun,x0,xt,y0,PointNum) %M

3、yEuler 用前向差分的欧拉方法解微分方程 %fun 表示f（x，y） %x0,xt表示自变量的初值和终值 %y0表示函数在x0处的值,其可以为向量形式 %PointNum表示自变量在x0,xt上取的点数 if nargin5 | PointNum dsolve(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,A,41,例题：用MATLAB的符号解法,求解常微分方程的特解：, dsolve(x*D

4、y+2*y-exp(x)=0,y(1)=2*exp(1),x) ans = (exp(x)*x-exp(x)+2*exp(1)/x2,A,42,例题：采用ODE45求解描述某刚性问题的微分方程：,function dy = 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-4 1e-4 1e-5); T,Y = ode45(rigid,0 12,0 1 1,

5、options); plot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.) legend(y1,y2,y3),A,43,A,44,例题：用MATLAB函数ode23,ode45,ode113,求解多阶常微分方程：,A,45,function dy=myfun03(x,y) dy=zeros(3,1) %初始化变量dy dy(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的三阶导数,% 用ode23

6、ode45 ode113解多阶微分方程 clear,clc x23,y23=ode23(myfun03,1,10,1 10 30); x45,y45=ode45(myfun03,1,10,1 10 30); x113,y113=ode113(myfun03,1,10,1 10 30); 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) %

7、以ode45为例作出函数以及其各阶导数图 legend(y,y一阶导数,y两阶导数) title(y,y一阶导数,y二阶导数函数图),A,46,A,47,例题：MATLAB在导热计算,传热过程涉及面广，数学模型复杂。计算过程中涉及到许多运算方法。导热虽是容易做数学处理的一种热量传递方式，但其过程往往涉及常微分、偏微分方程、线性（非线性）方程组的求解。,A,48,有一外径为4cm，内径为1.5cm，载有电流密度I为5000A/cm2的内冷钢质导体。导体单位时间发出的热量等于流体同时带走的热量，导体内壁面的温度为70。假定外壁面完全绝热。试确定，导体内部的温度分布；已知钢的导热系数k=0.38Kw

8、/mK，电导率=21011 km。,解：这是圆柱坐标中常物性一维稳态导热问题，结合本题具体条件导热微分方程式可简化为：,A,49,结合边界条件，可得这一导热问题的数学描述为：,此常微分方程的分析解，可调用MATLAB符号工具箱中的dsolve函数来实现。在命令窗口执行下面的代码： clear d_equat=D2t+Dt/r+131579=0; condition=（t0.0075）=70,D（t0.02）=0;%边界条件 t=dsolve（d_equat,condition,r） 程序执行结果：,A,50,程序执行结果：,即求出温度分布方程为:,A,51,工程上遇到的导热问题，往往由于物体的

9、几何形状复杂或边界条件难以描述，无法求出分析解，此时可采用数值方法进行求解。常微分方程（ODE）包括初值问题和边值问题两种，初值问题ODE的数值解法常调用ode45（）和ode23（）函数实现。 在MATLAB编辑器中编写函数BZ.m，存盘。 function BZ clear all clc a=0.0075;b=0.02;%r值的范围 solinit=bvpini（tlinspace（a,b,10）,70 80）；%用初始值对解初始化 sol=bvp4c（ODEfun,BCfun,solinit）；%解ODE方程 r=0.0075:0.00125:0.02 t=deva（lsol,r） function dtdr=ODEfun（r,t）%定义ODE方程函数 dtdr=（t2）;-131580-1/r*（t2）;