matlab 常微分方程数值解法

1、 物理学所研究的各种物质运动中，有许多物质运动的物理学所研究的各种物质运动中，有许多物质运动的过程是用常微分方程来描述的。过程是用常微分方程来描述的。例如例如, ,质点的加速运动，简谐振动等。质点的加速运动，简谐振动等。 dvFmdt22220d xxdt 简单问题可以求得解析解，多数实际问题靠数值求解。简单问题可以求得解析解，多数实际问题靠数值求解。 一阶常微分方程一阶常微分方程(ODE )初值问题初值问题 ：000( , ) ()ndyf x yxxxdxy xy()nnyy x 数值解法就是求数值解法就是求y(x)在某些分立的节点在某些分立的节点 xn 上的近似值上的近似值yn，用以近似

2、用以近似y(xn)ODE ：Ordinary Differential Equation2.1 简单欧拉简单欧拉(L.Euler, 1707-1783)方法。方法。00( , )()dyf y xdxy xy012nyyyy 欧拉数值算法就是由初值通过欧拉数值算法就是由初值通过递推求解递推求解，递推求解，递推求解就是从初值开始，后一个函数值由前一个函数值得到。关就是从初值开始，后一个函数值由前一个函数值得到。关键是构造递推公式。键是构造递推公式。2、欧拉近似方法、欧拉近似方法yyxxdd欧拉数值算法递推公式构造欧拉数值算法递推公式构造2.1.1 差分法差分法差分法就是用差商近似代替微商，即差分

3、法就是用差商近似代替微商，即代入微分方程得到：代入微分方程得到：()()( )( ( ),)( ( ), )yy xxy xxxy xf y xy x xxxxf对于等间隔节点对于等间隔节点11 n=0,1,2nnnnxxxhxxh 可以得到：可以得到：xn x0 x1 x2 . xn .y精确值精确值y(x0) y(x1) y(x2) . y(xn) .y近似值近似值y0 y1 y2 . yn .在在xn节点上，微分方程可以写为节点上，微分方程可以写为1()()() , xnnnny xy xf y xh作如下近似：作如下近似：( )nnyy t则得到欧拉解法递推公式的一般形式：则得到欧拉解

4、法递推公式的一般形式：1( , x )nnnnyyf yh具体求解过程为：具体求解过程为：1000( , x )yyf yh2111( , x )yyf yh3222( , x ) yyf yh简单欧拉方法程序简单欧拉方法程序function outx,outy=MyEuler(fun,x0,xt,y0,PointNum)%MyEuler 用前向差分的欧拉方法解微分方程用前向差分的欧拉方法解微分方程%fun 表示表示f（x，y）%x0,xt表示自变量的初值和终值表示自变量的初值和终值%y0表示函数在表示函数在x0处的值处的值,其可以为向量形式其可以为向量形式%PointNum表示自变量在表示自

5、变量在x0,xt上取的点数上取的点数if nargin5 | PointNum=0 %如果函数仅输入如果函数仅输入4个参数值，则个参数值，则PointNum默认值为默认值为100 PointNum=100;endif nargin 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例题：用例题：用MATLAB的符号解法的符号解法,求解常微分方程的特解：求解常微分方程的特解：1 20,2xxx

6、yyeye 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求解描述某刚性问题的微分方程：求解描述某刚性问题的微分方程：123121323123(0)0,(0)10.51(0)1yy yyyy yyyy yy 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 =

7、odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5);T,Y = ode45(rigid,0 12,0 1 1,options);plot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.)legend(y1,y2,y3)例题：用例题：用MATLAB函数函数ode23,ode45,ode113,求解多阶常微分方程求解多阶常微分方程：323232233(1)1,(1)10, (1)30,1,10 xd yd ydyxedxdxdxyyyx12232332333233xdyydxdyydxdyeyydxxxx1223333010000103230 x

8、ydYdyYdxdxyexxxfunction dy=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的三阶导数的三阶导数% 用用ode23 ode45 ode113解多阶微分方程解多阶微分方程clear,clcx23,y23=ode23(myfun03,1,10,1 10 30);x45,y45=

9、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) %以以ode45为例作出函数以及其各阶导数图为例作出函数以及其各阶导数图legend(y,y一阶导数一阶导数,y两阶导数

10、两阶导数)title(y,y一阶导数一阶导数,y二阶导数函数图二阶导数函数图)有一外径为有一外径为4cm，内径为，内径为1.5cm，载有电流密度，载有电流密度I为为5000A/cm2的内冷钢质导体。导体单位时间发出的热量等于流体同时带走的内冷钢质导体。导体单位时间发出的热量等于流体同时带走的热量，导体内壁面的温度为的热量，导体内壁面的温度为70。假定外壁面完全绝热。试。假定外壁面完全绝热。试确定，导体内部的温度分布；已知钢的导热系数确定，导体内部的温度分布；已知钢的导热系数k=0.38Kw/mK，电导率，电导率=21011 km。解：这是圆柱坐标中常物性一维稳态导热问题，结合本题具解：这是圆柱

11、坐标中常物性一维稳态导热问题，结合本题具体条件导热微分方程式可简化为：体条件导热微分方程式可简化为：结合边界条件，可得这一导热问题的数学描述为：结合边界条件，可得这一导热问题的数学描述为：此常微分方程的分析解，可调用此常微分方程的分析解，可调用MATLAB符号工具箱中的符号工具箱中的dsolve函数来实现。在命令窗口执行下面的代码：函数来实现。在命令窗口执行下面的代码：cleard_equat=D2t+Dt/r+131579=0;condition=（t0.0075）=70,D（t0.02）=0;%边界条件边界条件t=dsolve（d_equat,condition,r）程序执行结果：程序执行

12、结果：程序执行结果：程序执行结果：即求出温度分布方程为即求出温度分布方程为: :工程上遇到的导热问题，往往由于物体的几何形状复杂或边界条件难以描工程上遇到的导热问题，往往由于物体的几何形状复杂或边界条件难以描述，无法求出分析解，此时可采用数值方法进行求解。常微分方程（述，无法求出分析解，此时可采用数值方法进行求解。常微分方程（ODE）包括初值问题和边值问题两种，初值问题包括初值问题和边值问题两种，初值问题ODE的数值解法常调用的数值解法常调用ode45（）（）和和ode23（）函数实现。（）函数实现。在在MATLAB编辑器中编写函数编辑器中编写函数BZ.m，存盘。，存盘。function BZclear allclca=0.0075;b=0.02;%r值的范围值的范围solinit=bvpini（tlinspace（a,b,10）,70 80）；）；%用初始值对解初始化用初始值对解初始化sol=bvp4c（ODEfun,BCfun,solinit）；）；%解解ODE方程方程r=0.0075: