Matlab解微分方程ODEPDE

《Matlab解微分方程ODEPDE》由会员分享，可在线阅读，更多相关《Matlab解微分方程ODEPDE（43页珍藏版）》请在装配图网上搜索。

1、常微分方程:1 ODE解算器简介(ode*)2微分方程转换3刚性/非刚性问题(Stiff/Nonstiff)4隐式微分方程(IDE)5微分代数方程(DAE)6延迟微分方程(DDE)7边值问题(BVP)偏微分方程(PDEs)Matlab解法偏微分方程：1 一般偏微分方程组(PDEs)的命令行求解2特殊偏微分方程(PDEs)的PDEtool求解3陆君安偏微分方程的 MATLAB解法先来认识下常微分方程(ODE)初值问题解算器(solver)T,Y,TE,YE,IE = odesolver(odefun,tspan,yO,optionS sxint = deval(sol,xint)Matlab中提

2、供了以下解算器:解曽乐udEoivtr |趕44?* 5非m hrranztiffln：T- 4hi-怜式” Afta *ihiirif dr 2 NJ L. .tjLT K:=_耳.耳时T.岂x比ihiPtns.+ *.串K电十鬥亍卄韭用炜更曲拼狷拼 心iSJTAuu重麒lUJT展悝可RJT廿匕上”慎现 r+hrfii王古呼护離.占就茯用叫権常f. r九7-:*炖匕.确半囲寸馬馬H苕埼再+ 计/叭ME In i： 7-障佇ft ft+itfS住讯翳的聚鴉翠肝佥式，凶星冲田1*好什策 tW+ RClttHITmCft歯千言弹甘篇+ dtWltJRtttS悔些.r|t. 忌从上面的变换，我们注意

3、到，ODEs中所有因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式（比如上面的x（m）和y（n）不需要给它一个状态变量step3.根据上面选用的状态变量，写出所有状态变量的一阶微分的表达式X 賊+JS注意到，最高阶次的微分式的表达式直接就是stepl中的微分方程好，到此为止，一阶显式常微分方程组，变化顺利结束，接下来就是Matlab编程了。请记住上面的变化很重要，Matla中所有微分方程的求解都需要上面的变换。F面通过一个实例演示 ODEs的转换和求解己卸Apoll卫星的运动轨迹（兀刃满足下面的槪分方程组= 2y+ x-y = -2x-i-y讥并+旳乳4 = = 2兀2 + 兀3

4、口 JTj / q / 乜注意到，最高阶次的微分式的表达式直接就是stepl中的微分方程乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高阶次的微分式的表达式直接就是stepl中的微分方程其中，q = Jd + nF+yF = y/(x-y +/,出二打 82 45x(0) = 1.2,x(0) = 0JXO) = 0 j(0) = -l 04935751【解】 真是万幸，该 ODEs已经帮我们完成了 stepl，我们只需要完成 step2和step3了（1）选择一组状态变量乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高阶次的微分式的表达式直接就是s

5、tepl中的微分方程乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高阶次的微分式的表达式直接就是stepl中的微分方程（2）写出所有状态变量的一阶微分表达式乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高阶次的微分式的表达式直接就是stepl中的微分方程乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高阶次的微分式的表达式直接就是stepl中的微分方程2 = 2百+龙】一 “（可+“）r三”（兀一“）/勺乳4 = = 2兀2 + 兀3 口 JTj / q / 乜有了数学模型描述，则可以立即写出相应的Matlab代码了（这里我们优先

6、选择 ode45）1. x0=1.2;0;0;-1.04935751;%x0(i)对应与 xi的初值2. options=odeset(reltol,1e-8);% 该命令的另一种写法是options=odeset;options.reltol=1e-8;3. tic4. t,y=ode45(appollo,0,20,x0,options);%t 是时间点，y 的第 i 列对应 xi 的值，t 和 y 的行数相同5. toc6. plot(y(:,1),y(:,3)%绘制x1和x3,也就是x和y的图形7. title(Appollo卫星运动轨迹)8. xlabel( X)9. ylabel(Y

7、)10.10. function dx=appollo(t,x)11. mu=1/82.45;12. mustar=1-mu;13. r1= sqrt(x(1)+mu)A2+x(3)A2);14. r2=sqrt(x(1)-mustar)A2+x(3)A2);15. dx=x(2)16. 2*x(4)+x(1)-mustar*(x(1)+mu)/r1A3-mu*(x(1)-mustar)/r2A317. x(4)18. -2*x(2)+x(3)-mustar*x(3)/r1A3-mu*x(3)/r2A3;教程刚性/非刚性问题（Stiff/Nonstiff）的Matlab解法在工程实践中，我们经

8、常遇到一些ODEs其中某些解变换缓慢，另一些变化很快，且相差悬殊的微分方程，这就是所谓的刚性问题（Stiff），对于所有解的变化相当我们则称为非刚性问题（Nonstiff）。对于刚性问题一般不适合使用ode45这类函数求解。由于非刚性问题我们使用的多比较多，我们就不多说，下面主要讲解下Stiff看一个例子，考虑下面的微分方程” =0 04（1-孙）-（1-必）” + 0.0001（1-）2=-10000 + 30000-其中乃（0）二片3）二1 OJOO【解】 题目中已经帮我们完成了方程组转换，下面我们就直接编程了（1）编写微分方程函数1. odefun=(t,x)0.04*(1-x(1)-(

9、1-x(2)*x(1)+0.0001*(1-x(2)F22. -1e4*x(1)+3000*(1-x(2)F2;3. x0=0 1;4. tspan=0 100;太慢了)5. options=odeset(reltol,1e-6,abstol,1e-8);(2) 对于这个刚性问题，我们先使用ode45函数试试(反正我们是没有信心等待它,1. tic;t,y=ode45(odefun,tspan,x0,options);toc2. disp(ode45 计算的点数num2str(length(t) 换用ode15s函数试试看看1. tic;t2,y2=ode15s(odefun,tspan,x0

10、,options);toc2. disp(ode15s 计算的点数num2str(length(t2)(4)绘图看看结果1. figure(numbertitle,off,name,Stiff ODEs Demo by Matlabsky)2. subplot(121)3. title(Using ode45)4. plot(t,y)5. subplot(122)6. plot(t2,y2)7. title(Using ode15s)运行的结果如下，我们可以比较下ode45和ode15s的差距1. Elapsed time is 171.005688 seconds.2. ode45计算的点数

11、3569813. Elapsed time is 0.397287 seconds.4. ode15s计算的点数1881L,!Srng odsX&1U&mg心他甜-0?0B-C.B0.7-0.?2OS-0.50.S卩t-C403-03*0?-0.2010.1J-r.D(U1G08教程隐式微分方程(IDE)的Matlab解法上帝不会总是那么仁慈的，不是所有的ODEs都是可以直接显式的表达成下面的样子比如，下面的隐式微分方程组才+沙6ll+3x+zyl-z = 5那该如何办呢?在这里我们介绍三种解决方法，但是前两者是技巧方法，没有使用Mathworks专为IDE开发的ode15i函数：使用符号计算

12、函数 solve()求解出x和y的表达式，不就可以了吗？那好，我们下面试试看看解(1)求解表达式1.2. d2x,d2y=solve(d2x+2*dy*x-2*d2y,d2x*dy+3*dx*d2y+x*dy-y-5,d2x,d2y)%d2x 表示 x 的二阶导数，其它同理3.3. d2x =5.4. -2*(3*dy*x*dx-5+dy*x-y)/(3*dx+2*dy)7.5. d2y =9.6. (2*dyA2*x+5-dy*x+y)/(3*dx+2*dy)也就是说原来的ODEs可以重新写成亍=_2(3-5+7-7)3忑+2”y =也就是说stepl完成了，接下来的就比较简单了(2).选取

13、状态变量可二兀抵二X：碍二珂二卩(3) .写出状态变量的一阶微分表达式xj =严=_2張存1眄J+旺-可)(4) .写Matlab程序进行数值求解1. myode=(t,x)x(2)2. -2*(3*x(1)*x(2)*x(4)-5+x(4)-x(3)/(3*x(2)+2*x(4)3. x(4)4. (2*x(4)A2*x(1)+5-x(4)*x(1)+x(3)/(3*x(2)+2*x(4)5. t,y=ode45(myode,0,20,x0);但是solve()函数也不是万能的，对有些完全隐式，solve()可能压根没法求解出它们的具 体表达式，那此时该怎么办呢？天无绝人之路，车到山前必有路

14、，下面我们分析下数值解 法的本质。其实Matlab要求我们先将ODEs转化为显式的一阶微分方程组，就是为了便于计算状态变量一阶微分函数值(注意只需要它的值)！虽然现在我们没有解析的得到各个状态变量一阶 微分的表达式，但是我们只要求出每个点的对应的所有状态变量一阶微分值即可。这个思想就是在 Matlab没有提供ode15i()函数之前，大家唯一能使用的方法，下面我们以个例子说明下fsiii。)十 =-2 砂Ixx*y % cosfj;) = 3yx其中xO=XQO;l对于这样的ODEs难道你还想将它化为那个一阶显式微分方程组吗，别费劲了 ！【解】(1) 好，我们同样的需要一组状态变量西二兀 可二

15、瑞 兀二儿兀=yl(2) .写出状态变量的一阶微分。 只不过此时的x和y不再是显式的，我们使用fsolve进行数值求解花sin() + (x4 )3 = 一细码 + 码花 H珂；vj 筍 + cos(x4 *) =编写Matlab代码1. t,y=ode45(odefun,03,1 0 0 1);2. plot(t,y)3. Iegend(x,x,y,y”)4.4. function dy=odefun(t,x)5. fun=(dxy)dxy(1)*sin(x(4)+dxy(2F2+2*x(1)*x(3)-x(1)*dxy(1)*x(4)6. x(1)*dxy(1)*dxy(2)+cos(dx

16、y(2)-3*x(3)*x(2);7. options=optimset(display,off);y=fsolve(fun,x(1,3),options);%使用fsolve求解出x和 y10.dy=x(2);y(1);x(4);y(2);% 状态变量一阶微分值复制代码方法二也有它的缺点，每一个点的x“和y都需要独立计算， 假如说数据量很大的话，那个叫难熬呀，并且有时可能是由于方程或者参数配置问题fsolve还有会出现罢工现象！于是MathWorks为我们准备了 ode15i函数，不过这个 ode15i有点不好用，没有 ode45等 那么直接。ode15i()调用格式如下yOmod, yp0

17、mod,resnrm = decic(odefun,tO,yO,fixed_yO,ypO,fixed_ypO,options.)T,Y,TE,YE,IE = ode15i(odefu n,tspa n, y0mod,yp0mod,optio ns.)隐式ODEs不同于一般的显式 ODEs ode15i需要同时给出状态变量 及其一阶导数的初值（xO,dxO），它们不能任意赋值，只能有n个是独立的，其余的需要隐式方程求解（使用decic 函数），否则出现矛盾的初值条件。输入参数：odefun :微分方程，注意此时的函数变量有t（自变量），x（状态变量），dx（状态变量一阶导数）等三个tO :自变量

18、初值，必须与 tspan 相等y0 :状态变量初值，题目给出几个就写几个，没给的自己猜测一个fixed_yO :与yO的尺寸一样，表示对应位置的初值是给定的还是猜测的，如果给定则1，否则0ypO:状态变量一阶导数初值，其它同yOfixed_ypO :同理fixed_yO，不过此时对应的是 ypO 了options :其它选项，具体查看帮助tspan :自变量的范围，其中 tspan（1）必须与tO 一致输出参数：大部分同理ode45，这里就不具体说明了F面我们以实例说明，对于下面的IDE，求解x1口0) + (”呎=-2xy-xxny j+cosO-) = 3其中zO = l;ftO;l解（1

19、）选取状态变量西二兀 叼二瑞 為二儿旺=y(2)写出状态变量的一阶微分。只不过此时的x和y不再是显式的Xj 5111(X4)4-尸=-2矶可十兀/J 為珂花久*)=乡也也书写odefun，注意此时多了一个 dx变量，就是x的一阶导数1. odefun=(t,x,dx)dx(1)-x(2)2. dx(2)*si n(x(4)+dx(4F2+2*x(1)*x(3)-x(1)*dx(2)*x(4)3. dx(3)-x(4)4. x(1)*dx(2)*dx(4)+cos(dx(4)-3*x(3)*x(2);(4)求解dx01. t0=0;%自变量初值2. %对于x0和dx0中的题目给出的初值，如实写，

20、没有给出的任意写3. x0=1 0 0 1;%本题初值x0的都给出了，所以必须是这个4. %初值x0中题目给出的，对应位置填1,否则为05. fix_x0=ones(4,1);%本题中x0都给出了，故全为 16. dx0=0 0 1 1;%本题中初值dx0 一个都没有给出，那么全部任意写7. %初值dx0中题目给出的，对应位置填 1，否则为08. fix_dx0=zeros(4,1);%本题中dx0 一个没有给出，故全部为 09. x02,dx02=decic(odefun,t0,x0,fix_x0,dx0,fix_dx0);(5) 求解微分方程1. solution=ode15i(odefu

21、n,0 2,x02,dx02)%x02 和 dx02 是由 decic 输出参数虽然我们上面的分析是完全正确，但是好像Matlab就是罢工，郁闷死了其实对一些简单的 ODEs压根没有必要通过 decic函数来求解未知的初值，我们可以直接人工算出x02和dx02，根据下面的式子花 011（驻）+ E 忙=一2誥内 + X2 *X4我们演示下，根据题目x0=x1,x2,x3,x4=1 0 0 1，由第一个式子可以得出x1 =0,第三个式子有x3 =1，再联立第二四两个式子后有x2 =0，x4 =1，故dx0=x1 ,x2 ,x3 ,x4 =0 0 1 1看了下，隐式 ODEs数值求解的确麻烦的些，

22、但是和前面的ode45没有本质的区别，decic函数就是为了求解那些隐式中的状态变量的一阶导数，与我们方法二中的很相似，只不过 这里直接调用罢了。其实不管隐式还是显式 ODEs甚至DAE都是可以使用 ode15i求解，只不过此时将显式当 隐式处理，或者说将代数约束看成一个隐式罢了！这一条等到我们在介绍 DAE的时候再验证！当然大部分人都不会傻到这么做，有现成的函数不用，何必自找麻烦呢！但是对于了解和 学习，我们可以试试。教程微分代数方程(DAE)的Matlab解法所谓微分代数方程，是指在微分方程中，某些变量满足某些代数方程的约束。假设微分方 程的更一般形式可以写成前面所介绍的ODEs数值解法主

23、要针对能够转换为一阶常微分方程组的类型，故DAE就无法使用前面介绍的常微分方程解法直接求解，必须借助DAE的特殊解法。其实对于我们使用 Matlab求解DAE时，却没有太大的改变 只需增加一个 Mass参数即可。 描述f(t,x)的方法和普通微分方程完全一致。、，、，、厂注意： ode15i 没法设置 Mass 属性，换句话说除了 ode15i 外其他 ode 计算器都可以求解DAEs问题1. 当 M(t,y) 非奇异的时候， 我们可以将微分方程等效的转换为 y=inv(M)*f(t,y) ，此时就是一个普通的ODE当然我们可以将它当成 DAEs处理)，对任意一个给定的初值条件都有唯一的解2.

24、 当m(t,y)奇异时，我们叫它为 DAEs(微分代数方程)，DAEs问题只有在同时提供状态变量初值y0和状态变量一阶导数初值pyO，且满足M(tO,yO)*ypO=f(tO,yO)时才有唯一解，假如不满足上面的方程，DAEs解算器会将提供的yO和pyO作为猜测初始值，并重新计算与提供初值最近的封闭初值3. 质量矩阵可是一个常数矩阵 (稀疏矩阵 )，也可以是一个自定义函数的输出。但是 ode23s只能求解Mass是常数的DAEs4. 对于Mass奇异的DAEs问题，特别是设置 MassSingular为yes时，只能ode15s和ode23t 解算器5. 对于DAE我们还有几个参数需要介绍a.

25、 Mass :质量矩阵，不说了，这个是DAE的关键，后面看例子就明白了b. MStateDependence ：质量矩阵 M(t,y) 是否是 y 的函数，可以选择 none|weak|strong ， none 表示 M与 y 无关，weak禾口 strong 者E表示与 y 相关c. MvPattern :注意这个必须是稀疏矩阵， S(i,j)=1 表示 M(t,y) 的第 i 行中任意元素都与 第 j 个状态变量 yi 有关，否则为 Od. MassSingular :设置 Mass矩阵是否奇异，当设置为yes时，只能使用 ode15s和ode23te. InitialSlope:状态变

26、量的一阶导数初值yp0，和y0具有相同的size，当使用ode15s和 ode23t 时，该属性默认为 O面我们以实例说明，看下面的例子，求解该方程的数值解珂 =0-2s)十 叼珂 十0 2珂屯 = 2屯;t - 5x32r? - 2才再+吃+旳=1其中匝() 0 & 旳()=忑(0) = 0.1解真是万幸，选取状态变量和求状态变量的一阶导数等，微分方程转换工作，题目已经帮我 们完成。DAE可是细心的网友会发现，最后一个方程不是微分方程而是一个代数方程(这就是为什叫的原因)，其实我们可以将它视为对三个状态变量的约束。(1) 用矩阵形式表示出该DAEs1 0 CT -0 2阳+ %再+ 0.引再

27、0 1 02珂X： -5%； -0 0 0%! + 2 4-1编写Matlab代码1. odefun=(t,x)-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);2. 2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);3. x(1)+x(2)+x(3)-1;% 微分方程组4. M=1 0 0;0 1 0;0 0 0;% 质量矩阵5. options=odeset(mass,M);%对以 DAE 问题，mass属性必须设置6. x0=0.8;0.1;0.1;% 初值7. t,x=ode15s(odefun,0 20,x0,options);% 这里好像不能

28、使用 ode458. figure(numbertitle,off,name,DAE demo by Matlabsky)plot(t,x)10.Iegend(x1(t),x2(t),x3(t)复制代码在介绍隐式ODEs时，我们说了 ode15i也可以求解DAEs原理就是将约束看成一个隐式 好，下面我们看看，验证下所说的：1. odefun=(t,x,dx)dx(1)+0.2*x(1)-x(2)*x(3)-0.3*x(1)*x(2);2. dx(2)-2*x(1)*x(2)+5*x(2)*x(3)+2*x(2)A2;3. x(1)+x(2)+x(3)-1;4. %状态变量初值，题目中给出5.

29、x0=0.8;0.1;0.1;6. %该初值全部给定，按理说应该全部为1，但是记住方程中有一个约束，故其实只有两个独立初值7. x0_fix=1;0;1;%随意一个改为 0都可以，比如0;1;1或者1;1;08. %状态变量一阶微分初值，题目没有给岀，可以任意写dx0=1;1;1;9. %该初值一个都没有给出，故全部为010. dxO_fix=O;O;O;11. %时间变量的初值12. tO=O;13. %计算相容初值14. xO,dxO=decic(odefun,O,xO,xO_fix,dxO,dxO_fix);15. xO,dxO % 相容初始条件16. res=ode15i(odefun

30、,O,2O,xO,dxO);17. figure(numbertitle,off,name,DAE demo by Matlabsky)18. plot(res.x,res.y)其实，有些简单DAEs是可以转换为ODEs的，由于代数方程只是状态变量的一个约束， 假如约束充分的话，我们就可以使用消参的思想，将约束反代回ODEs中，这样就将那个约束方程消去了，最后只剩下那些微分方程了。对于本例，我们只要通过第三个方程表示出 再将x3反代回前两个微分方程中则有 = -0.2 +Xj-兀2)+ 0. 3阿兀2花二 2 珂电5a(1 Aj a) 2aj这个就是一个正常的ODEs了，好，下面我们就重新使用

31、Matlab求解下看:1. x0=0.8;0.1;2. odefun=(t,x)-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2)3. 2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2);4. t,x=ode45(odefun,0,20,x0);5. figure(numbertitle,off,name,DAE demo by Matlabsky)6. plot(t,x)7. Iegend(x1(t),x2(t)教程延迟微分方程(DDE)的Matlab解法延迟微分方程的一般形式为：兀(f) = 菇xQ-巧)/(一%),盂一耳

32、)其中taui 0为状态变量x(t)的延迟常数。Matlab中提供了求解这类方程的隐式Runge-Katta算法dde23()函数，可以直接求解延迟微分方程，该函数的格式为：sol = dde23(ddefun,tau,history,tspan,options )输入参数：ddefun :描述延迟微分方程的句柄tau=tau1,taou, .taun:history :描述t (/)(2) 计算每个状态变量的一阶微分珂(Z) = 1 3x(f)疋2 ( 1 0.2 f ( 0.5) Xj 0.5) 宅可O =码)花(f)二 4 吗() 一 (f) 一 3 也()编写延迟函数。从第一个方程可以

33、看出x2延迟1, x1延迟0.5，我们取tau=1 0.51. function dx=ddefun(t,x,z)2. %z(i,j)表示状态变量xj延迟了tau(i)3. tau仁z(:,1);%第一列表示延迟了tau(1)=1的所有状态变量4. tau2=z(:,2);5. %x2延迟了 1,故表示为z(2,1)或者tau1(2)6. %x1延迟了 0.5，故表示为z(1,2)或者tau2(1)7. dx=1-3*x(1)-tau1(2)-0.2*tau2(1F3-tau2(1)8. x(3)9. 4*x(1)-2*x(2)-3*x(3);(4)编写主调函数1. tau=1 0.5;2.

34、history=0 0 0;3. tspan=0 10;4. sol=dde23(ddefun,tau,history,tspan);5. figure(numbertitle,off,name,DDE demo by Matlabsky)6. plot(sol.x,sol.y(2,:)教程边值问题(BVP)的 Matlab解法关于边值微分方程问题相关连接有：【非线性微分方程边值问题打靶算法】参见【线性微分方程边值问题打靶算法】参见 【常微分方程 - 边值问题六阶精度解算器 (BVP6C)】 这里有好多 bvp4c 使用例子，不过 是英文的 到现在为止，前面讲到的都是初值问题，也就是说 Mat

35、lab 的 ode* 系列解算器，默认将 tspan(1) 作为初值条件时的 t ，比如你将初值条件换为 x(2)=x (2)=0 ， 那么 tspan(1) 就必 须是 2但是工程应用中我们经常遇到边值问题， 这些是那些 ode* 函数无能为力的， 当然我们可以 自己编写函数求解 ( 比如 shooting) ，但是那个毕竟不是某些人能力所及的， 还好 Matlab 中 提供了 bvp 解算器。solinit = bvpinit(x, yinit,params)sol = bvpsolver(odefun,bcfun,solinit,options )由于边值问题可能有多解， 为了便于我们确

36、定那个解是我们需要的， 所以必须使用 bvpinit 函数对初值进行估计解算器 (bvpsolver) ： Matlab 中提供了 bvp4c 和 bvp5c ，后者误差控制更好些输入参数：x ：需要计算的网格点，相当于 ode* 的 tspanyinit ：状态变量的初始猜测的值，可以是具体值，也可以是函数，类似与 ode* 的 x0params:微分方程中的其它未知参数的初始猜测值odefun ：描述边值问题微分方程的函数句柄bcfun :边值函数，一般是双边值 (x 的上下限即认为两个边界 ) ，但也支持多边值 (具体看帮 助)solinit :由 bvpinit 生成的初始化网格opt

37、ions : BVP解算器优化参数，可以通过 bvpset设置，具体参数查看帮助输出参数：大部分同理ode45，这里不详细介绍了入也是未好下面我们以例子说事，考察包含未知参数入的二阶微分方程，求解入，由于知的(这就是在params中说到的其他位置参数)，故需要 三个边值条件y+ (2 - 2q cos 2x) = 0 其中q = 5,y0) =y(7f) = 0, y (0) = 1解(1)选用状态变量(2)求每个状态变量的一阶微分 = y= x3花 1 =丸 + 2q cos 2： j(3) 编程实现求解1. q=5;2. lambda = 15;%未知参数猜测值3. x=linspace(

38、0,pi,10);% 需要计算的点4. yinit=(x)cos(4*x);-4*sin(4*x);% 使用函数估计目标值5. solinit = bvpinit(x,yinit,lambda);% 生成计算网格6. odefun=(x,y,lambda)y(2);-(lambda - 2*q*cos(2*x)*y(1);% 边值微分方程7. bcfun=(ya,yb,lamdba)ya(2);yb(2);ya(1)-1;% 边界条件8. sol = bvp4c(odefun,bcfun,solinit);9. fprintf(The fourth eigenvalue is approxim

39、ately %7.3f.n,sol.parameters)10. xint = linspace(0,pi);11. Sxint = deval(sol,xint);% 计算 xint 对应的 y 值12. figure(name,BVP Demo by Matlabsky, numbertitle , off )13. plot(xint,Sxint(1,:)% 绘制 xint 和 x1=y 的图形，即 x 与 y14. axis(0 pi -1 1.1)15. title(Eigenfunction of Mathieus equation.)16. xlabel(x)17. ylabel

40、(solution y)已;un lion of弓“川前沖般偏微分方程组(PDEs)的MATLAB求解一、pdepe()函数说明MATLAB语言提供了 pdepe()函数，可以直接求解一般偏微分方程(组)，它的调用格式为1. sol=pdepe(m,pdefun,pdeic,pdebc,x,t)【输入参数】pdefun：是PDE的问题描述函数，它必须换成下面的标准形式0 （和浮）申=计（兀覚见舟）十卫,申）（式1）办 atoxaxdx这样，PDE就可以编写下面的入口函数1. c,f,s=pdefun（x,t,u,du）c,f,s这三个函m,x,t就是对应于（式1）中相关参数，du是u的一阶导数

41、，由给定的输入变量即可表示出出 数pdebc：是PDE的边界条件描述函数，必须先化为下面的形式戸（心3 +0（和忌）护r） = 0于是边值条件可以编写下面函数描述为1. pa,qa,pb,qb=pdebc（x,t,u,du）其中a表示下边界，b表示下边界pdeic：是PDE的初值条件，必须化为下面的形式股我们使用下面的简单的函数来描述为1. uO=pdeic(x)m,x,t:就是对应于(式1)中相关参数【输出参数】sol：是一个三维数组，sol(:,:,i)表示ui的解，换句话说 uk对应x(i)和t(j)时的解为sol(i,j,k)通过sol，我们可以使用pdeval()直接计算某个点的函数

42、值1、PDEs问题实例讲解试求解下面的偏微分鲁024勢-叫-吗) 竽17券+盹宀)其中,F(x)=严-11.46r且满足初始条件 :-，及边界条件学3丿)=眄(6Q = o.坷(1丿)=h学1丿)=oOXOX【解】（1 ）对照给出的偏微分方程，根据标注形式，则原方程可以改写为（旳一旳） F （叭一妁）注意：1;1那里是“.*而不是“，”表示对应元素相乘，而不是矩阵的相乘，不要问为什么！这就是Matlab规定的必须转化的标准形式可见m=0 ,且1. %目标PDE函数2. function c,f,s=pdefun (x,t,u,du)3. c=1;1;4. f=0.024*du(1);0.17*

43、du(2);5. temp=u(1)-u(2);6. s=-1;1.*(exp(5.73*temp)-exp(-11.46*temp);(2 )边界条件改写为下边界o上边界1. %边界条件函数2. function pa,qa,pb,qb=pdebc(xa,ua,xb,ub,t)3. %a表示下边界，b表示上边界4. pa=0;ua(2);5. qa=1;0;6. pb=ub(1)-1;0;7. qb=0;1;(3 )初值条件改写为1. %初值条件函数2. function u0=pdeic(x)3. u0=1;0;(4)最后编写主调函数1. clc2. x=0:0.05:1;3. t=0:0

44、.05:2;4. m=0;5. sol=pdepe(m,pdefun,pdeic,pdebc,x,t);by Matlabsky)6. figure(numbertitle,off,name, P DE Demo 7. subplot(211)8. surf(x,t,sol(:,:,1)9. title(The Solution of u_1)10. xlabel(X)11. ylabel(T)12. zlabel(U)13. subplot(212)14. surf(x,t,sol(:,:,2)15. title(The Solution of u_2)16. xlabel(X)ylabel

45、(T)18. zlabel(U) All典型偏微分方程的描述下面我们先了解下三个典型的二阶PDE，然后介绍pdetool，至于命令行我们就免了，真的很累人，如果的确需要的话，那就让 Matalb直接生成就好了。(1) 椭圆型 偏微分方程的一般形式为div(cVu) +a*u = 即弋(+?其中3 J为给定的函数或者常数(2) 抛物线型偏微分方程的一般形式护旦+十厶)十几况水字一出卩2讷+口 即 厂史(三+ ,dt迦/其中乩G 口必须是常数(3) 双曲线型 偏微分方程的一般形式%d v一div(cV) -au =即2 .-T|2幷2d牛一-(： *( 4三+如+匕*盘二丁(工丄)审 対禺/东7八

46、其中Rc/必须是常数一应卩（/7盘）4卫也J =刃卅占冷肚即从上面可以看出， 三类典型二阶偏微分方程的区别在于u对t的导数阶次。椭圆型PDEs中,和f可以是给定的函数或者常数，但是其它两类必须都是常数。c、a、dMATLAB是采用有限元的方法求解各种PDE。MATLAB为我们提供一个 pdetool的交互界面,解二元偏微分 u（x1,x2）（注意只能求解二元）。方程的参数由a、c、d和f确定，求解域由图形确定，求解域确定好后，需要对求解域进行栅格化 个是自动，我只要进行相关设置进行网格控制就可以了）。可以求（这5#dpdetool可以执行解的动态显示，也就是根据时间参数t，绘制pde问题的求解过程y动画播放具体设置如下：m【Plot】-【Parameters选中【Animation】，点击后面的【Options】，设置播放 速度和次数，比如6fps表示每秒6帧c【Plot-【Export Movie 输入动画保存的变量名，比如M(4) 在 Command Windows 中直接输入 movie(M) 即可播放使用movie2ve(M, demo.avi命令可以将动画保存为avi文件