Matlab解微分方程(ODE+PDE)

常微分方程：1ODE解算器简介(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,y0,options)sxint=deval(sol,xint)Matlab中提供了以下解算器：输入参数：odefun：微分方程的Matlab语言描述函数，必须是函数句柄或者字符串，必须写成Matlab规范格式(也就是一阶显示微分方程组)，这个具体在后面讲解tspan=[t0tf]或者[t0,t1,…tf]：微分变量的范围，两者都是根据t0和tf的值自动选择步长，只是前者返回所有计算点的微分值，而后者只返回指定的点的微分值，一定要注意对于后者tspan必须严格单调，还有就是两者数据存储时使用的内存不同(明显前者多)，其它没有任何本质的区别y0=[y(0),y’(0),y’’(0)…]：微分方程初值，依次输入所有状态变量的初值，什么是状态变量在后面有介绍options：微分优化参数，是一个结构体，使用odeset可以设置具体参数，详细内容查看帮助输出参数：T：时间列向量，也就是ode**计算微分方程的值的点Y：二维数组，第i列表示第i个状态变量的值，行数与T一致在求解ODE时，我们还会用到deval()函数，deval的作用就是通过结构体solution计算t对应x值，和polyval之类的很相似！参数格式如下：sol：就是上次调用ode**函数得道的结构体解xint：需要计算的点，可以是标量或者向量，但是必须在tspan范围内该函数的好处就是如果我想知道t=t0时的y值，不需要重新使用ode计算，而直接使用上次计算的得道solution就可以[教程]微分方程转换为一阶显示微分方程组方法好，上面我们把Matlab中的常微分方程(ODE)的解算器讲解的差不多了，下面我们就具体开始介绍如何使用上面的知识吧！现实总是残酷的，要得到就必须先付出，不可能所有的ODE一拿来就可以直接使用，因此，在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步，借助状态变量将微分方程组化成Matlab可接受的标准形式(一阶显示常微分方程)！如果ODEs由一个或多个高阶微分方程给出，则我们应先将其变换成一阶显式常微分方程组！下面我们以两个高阶微分方程构成的ODEs为例介绍如何将之变换成一个一阶显式常微分方程组。step1.将微分方程的最高阶变量移到等式的左边，其他移到右边，并按阶次从低到高排列，假如说两个高阶微分方程最后能够显式的表达成如下所示：我们说过现实总是残酷的，有时方程偏偏是隐式的，没法写成上面的样子，不用担心Matlab早就为我们想到了，这个留在后面的隐式微分方程数值求解中再详细讲解！step2.为每一阶微分式选择状态变量，最高阶除外从上面的变换，我们注意到，ODEs中所有因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式(比如上面的x(m)和y(n))不需要给它一个状态变量step3.根据上面选用的状态变量，写出所有状态变量的一阶微分的表达式注意到，最高阶次的微分式的表达式直接就是step1中的微分方程好，到此为止，一阶显式常微分方程组，变化顺利结束，接下来就是Matlab编程了。请记住上面的变化很重要，Matla中所有微分方程的求解都需要上面的变换。下面通过一个实例演示ODEs的转换和求解【解】真是万幸，该ODEs已经帮我们完成了step1，我们只需要完成step2和step3了(1)选择一组状态变量(2)写出所有状态变量的一阶微分表达式(4)有了数学模型描述，则可以立即写出相应的Matlab代码了(这里我们优先选择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')10.11.functiondx=appollo(t,x)12.mu=1/82.45;13.mustar=1-mu;14.r1=sqrt((x(1)+mu)^2+x(3)^2);15.r2=sqrt((x(1)-mustar)^2+x(3)^2);16.dx=[x(2)17.2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^318.x(4)19.-2*x(2)+x(3)-mustar*x(3)/r1^3-mu*x(3)/r2^3];[教程]刚性/非刚性问题(Stiff/Nonstiff)的Matlab解法在工程实践中，我们经常遇到一些ODEs，其中某些解变换缓慢，另一些变化很快，且相差悬殊的微分方程，这就是所谓的刚性问题(Stiff)，对于所有解的变化相当我们则称为非刚性问题(Nonstiff)。对于刚性问题一般不适合使用ode45这类函数求解。由于非刚性问题我们使用的多比较多，我们就不多说，下面主要讲解下Stiff看一个例子，考虑下面的微分方程【解】题目中已经帮我们完成了方程组转换，下面我们就直接编程了(1)编写微分方程函数1.odefun=@(t,x)[0.04*(1-x(1))-(1-x(2))*x(1)+0.0001*(1-x(2))^22.-1e4*x(1)+3000*(1-x(2))^2];3.x0=[01];4.tspan=[0100];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))])(3)换用ode15s函数试试看看1.tic;[t2,y2]=ode15s(odefun,tspan,x0,options);toc2.disp(['ode15s计算的点数'num2str(length(t2))])(4)绘图看看结果1.figure('numbertitle','off','name','StiffODEsDemo—byMatlabsky')2.subplot(121)3.title('Usingode45')4.plot(t,y)5.subplot(122)6.plot(t2,y2)7.title('Usingode15s')运行的结果如下，我们可以比较下ode45和ode15s的差距1.Elapsedtimeis171.005688seconds.2.ode45计算的点数3569813.Elapsedtimeis0.397287seconds.4.ode15s计算的点数188[教程]隐式微分方程(IDE)的Matlab解法上帝不会总是那么仁慈的，不是所有的ODEs都是可以直接显式的表达成下面的样子比如，下面的隐式微分方程组那该如何办呢？在这里我们介绍三种解决方法，但是前两者是技巧方法，没有使用Mathworks专为IDE开发的ode15i函数：================方法一================使用符号计算函数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.4.d2x=5.6.-2*(3*dy*x*dx-5+dy*x-y)/(3*dx+2*dy)7.8.d2y=9.10.(2*dy^2*x+5-dy*x+y)/(3*dx+2*dy)也就是说原来的ODEs可以重新写成也就是说step1完成了，接下来的就比较简单了(2).选取状态变量(3).写出状态变量的一阶微分表达式(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)^2*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()可能压根没法求解出它们的具体表达式，那此时该怎么办呢？天无绝人之路，车到山前必有路，下面我们分析下数值解法的本质。其实Matlab要求我们先将ODEs转化为显式的一阶微分方程组，就是为了便于计算状态变量一阶微分函数值(注意只需要它的值)！虽然现在我们没有解析的得到各个状态变量一阶微分的表达式，但是我们只要求出每个点的对应的所有状态变量一阶微分值即可。这个思想就是在Matlab没有提供ode15i()函数之前，大家唯一能使用的方法，下面我们以一个例子说明下对于这样的ODEs难道你还想将它化为那个一阶显式微分方程组吗，别费劲了！【解】(1)好，我们同样的需要一组状态变量(2).写出状态变量的一阶微分。只不过此时的x’’和y’’不再是显式的，我们使用fsolve进行数值求解(3)编写Matlab代码1.[t,y]=ode45(@odefun,[03],[1001]');2.plot(t,y)3.legend('x','x''','y','y''')4.5.functiondy=odefun(t,x)6.fun=@(dxy)[dxy(1)*sin(x(4))+dxy(2)^2+2*x(1)*x(3)-x(1)*dxy(1)*x(4)7.x(1)*dxy(1)*dxy(2)+cos(dxy(2))-3*x(3)*x(2)];8.options=optimset('display','off');9.y=fsolve(fun,x([1,3]),options);%使用fsolve求解出x''和y''10.dy=[x(2);y(1);x(4);y(2)];%状态变量一阶微分值复制代码================方法三================方法二也有它的缺点，每一个点的x''和y''都需要独立计算，假如说数据量很大的话，那个叫难熬呀，并且有时可能是由于方程或者参数配置问题fsolve还有会出现罢工现象！于是MathWorks为我们准备了ode15i函数，不过这个ode15i有点不好用，没有ode45等那么直接。ode15i()调用格式如下[y0mod,yp0mod,resnrm]=decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options...)[T,Y,TE,YE,IE]=ode15i(odefun,tspan,y0mod,yp0mod,options...)隐式ODEs不同于一般的显式ODEs，ode15i需要同时给出状态变量及其一阶导数的初值(x0,dx0)，它们不能任意赋值，只能有n个是独立的，