如何使用MATLAB求解微分方程(组).ppt

TMU_BME_2013,Topic: 如何使用MATLAB求解常微分方程（组）,a.What ?,微分方程指描述未知函数的导数与自变量之间的关系的方程。未知函数是一元函数的微分方程称作常微分方程。未知函数是多元函数的微分方程称作偏微分方程。 MATLAB（matrix&laboratory）意为矩阵工厂（矩阵实验室）.MATLAB是美国MathWorks公司出品的商业数学软件，提供高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。,Where ?,ODE,工程控制,航空航天,金融分析,医学生理,Where ?,MAT LAB,算法开发,数据可视化,数据分析,数值计算,When ?,当对问题进行建模后，有常微分方程需要求解时。 在生物建模中，经常需要求解常微分方程。如药物动力学的房室模型的建模仿真。,How ?,数值解,在数学中，刚性方程是指一个微分方程，其数值分析的解只有在时间间隔很小时才会稳定，只要时间间隔略大，其解就会不稳定。 目前很难去精确地去定义哪些微分方程是刚性方程，但是大体的想法是：这个方程的解包含有快速变化的部分。,数值解是采用如有限元方法、 数值逼近方法、插值方法等计算方法得到的解。只能利用数值计算的结果， 而不能随意给出自变量并求出计算值。 当无法藉由微积分技巧求得解析解时，这时便只能利用数值分析的方式来求得其数值解了。实际情况下，常微分方程往往只能求解出其数值解。,数值解？刚性方程？,如何调用？,y=dsolve(e1,e2,.,c1,c2,.,v) 其中e1,e2,.为微分方程或微分方程组； c1,c2,.,是初始条件或边界条件； v是独立变量，默认的独立变量是t； y 返回解析解。如果没有初始条件，则求出通解，如果有初始条件，则求出特解。 用字符串表示常微分方程，自变量缺省时为t,导数用D表示微分。y的2阶导数用D2y表示，依此类推。,如何调用？,T,Y,TE,YE,IE=solver(odefun,tspan,y0,options) 其中solver为ode23、ode45、ode113、ode15s、ode23s、ode23t、ode23tb 函数； odefun 是函数句柄； tspan 微分定义区间； y0 为初值行矩阵； T 值是t序列（为列向量）； Y 值是微分方程的解Y在各点t的值（为列向量）； TE 表示事件发生时间，可缺省； YE 表示事件解决时间，可缺省； IE 表示事件消失时间，可缺省； options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等，可缺省。,使用ODE?时如何编写微分方程 ？,方式一：带额外参数，使用时需对参数进行赋值 function odefun(t,x,flag,R,L,C) %用flag说明R、L、C为变量 xdot=zeros(2,1); %表明其为列向量 xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t); xdot(2)=1/C*x(1); end,方式二：无额外参数 function dC=odefun(t,C) dC= 0.1*C(3)+37.5-2.52*C(1); %直接书写列矩阵 1.68*C(1)-0.01*C(2); 0.004*C(2)-0.1*C(3); end,Examples,E.g.1 求下列微分方程组的通解 x=2x-3y+3z y=4x-5y+3z z=4x-4y+2z,输入命令： x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z, Dz=4*x4*y+2*z,t) x=simple(x); y=simple(y); z=simple(z); %化简结果,运行结果为： x=(c1-c2+c3+c2e -3t-c3e-3t)e2t y= - c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z=(-c1 e-4t+c2e-4t+c1-c2+c3)e2t,Examples,E.g.2 求解y= - t*y + et*y +3sin(2t)。 已知3.9t4.0，y(0)=2，y(0)=8,函数文件odefun.m的建立 function y=odefun(t,x) y=zeros(2,1); %列向量 y(1)=x(2); y(2)= -t*x(1)+exp(t)*x(2)+3*sin(2*t); end 求解程序关键步骤 t,y=ode45(odefun,3.9 4.0,2 8),程序执行结果,E.g.3 求解人体中不同形式的碘浓度的三房室模型。已知碘在三房室之间的转换速率为：k21=0.84/d，k01=1.68/d，k32=0.01/d，k13=0.08/d，k03=0.02/d，f10=150g/d，f30=0g/d。初始条件为：x1(0)=81.2g，x2(0)=6821g，x3(0)=682g，v1=4L, v2=2L, v3=5L。仿真时间为30d。,Examples,通过列写方程组，建立odefun.m文件即方程组文件 function dC=odefun(t,C) dC=0.1*C(3)+37.5-2.52*C(1); 1.68*C(1)-0.01*C(2); 0.004*C(2)-0.1*C(3); end 程序关键部分列写如下： t,C=ode15s(fun1,0,30,20.3,3410.5,136.4),E.g.3结果,Examples,E.g.4 求解方程y+1000(y2-1)y+y=0。已知