MATLAB实验4-微分方程求解.ppt

实验名称：微分方程求解 实验目的：熟悉并掌握常微分方程数值求解的基本原理、方法；掌握常微分方程模型的MATLAB解析解法与数值解法。 实验任务： 1、P139，习题6 2、P140，习题13. 实验4 附录一：基于MATLAB的微分方程求解 一、微分方程的解析运算 y=diff(fun,x) %函数fun对x求导数 1、函数求导的解析运算 y=diff(fun,x,n) %函数fun对x求n阶导数 例： 求 syms x f=sin(x)/(x^2+4*x+3);f1=diff(f,x,4);pretty(f1) 2、微分方程的解析运算 MATLAB符号运算工具箱提供了一个线性常系数微分方程的求解函数dsolve()，该函数允许用字符串的形式描述微分方程及初值、边值条件，最终得出微分方程的解析解。调用格式： y=dsolve(f1,f2, …,fm,’t’) 其中fi既可以描述微分方程，又可以描述初始条件或边界条件。 用D4y 表示，已知条件 ,可用 D2y(2)=3表示。 例1： ，求下列微分方程的通解 syms t y u=exp(-5*t)*cos(2*t+1)+5; uu=diff(u,t,2)+4*diff(u,t)+2*u; dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=uu) 若考虑初始条件： 则： y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y=uu,y(0)=3,Dy(0)=2,D2y(0)=0,D3y(0)=0) 例2：求下列线性微分方程组 [x,y]=dsolve(D2x+2*Dx=x+2*y-exp(-t),Dy=4*x+3*y+4*exp(-t)) 二、微分方程数值求解 1、Euler法，梯形公式，预估-校正法，线性多步法，Runge-Kutta等。 2、数值求解的MATLAB函数 (1)、ode45()，是非刚性问题大部分场合的首选算法。是单步算法，采用4，5阶Runge-Kutta方程，也称为RKF格式。 (2)、ode23()，适用于非刚性问题精度较低的情形，是单步算法，采用2，3阶Runge-Kutta方程。 (3)、ode113()，适用于非刚性问题，采用Adams多步法。 (4)、ode23t()，适用于适度刚性问题，采用梯形算法。 (5)、ode15s()，适用于刚性问题，若ode45失败时可尝试使用，采用多步法，Gear’s反向数值微分，精度中等。 (6)、ode23s()，适用于刚性问题，采用单步法，2阶Rosebrock算法，精度低。 下面，以ode45为例说明函数的用法。调用格式： [t,x]=ode45(Fun,ts,x0) [t,x]=ode45(Fun,ts,x0,options) 其中，Fun是由微分方程或方程组写成的M文件，ts描述方程的求解区间，如果是二维向量[t0,tf]表示自变量初值t0和终值tf，如果是高维向量[t0,t1,…,tf](或者t0:k:tf),则表示输出这些节点的值。x0为函数的初值，options为控制选项，可以通过odeset()函数获取。 常用的控制参数为： RelTol~相对误差容许上限，默认值为0.001，为了保证较高的精度可以减少该值 。 AbsTol~为一个向量，其分量表示每个状态变量允许的绝对误差，其默认值为10^(-6)，可以改变求解精度。 例如，相对误差要求为10^(-7)时 options=odeset(‘RelTol’,1e-7); 例：设著名的Lorenz模型的状态方程表示为 其中，设 其初始值为 ，求解该微分方程组。 function xdot=lorenzeq(t,x) xdot=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)]; clear;x0=[0;0;1e-10]; [t,x]=ode45(@lorenzeq,[0,100],x0); subplot(1,2,1) plot(t,x); title(‘状态变量的时间响应图); subplot(1,2,2) plot3(x(:,1),x(:,2),x(:,3));title(‘相空间三维图); axis([10 42 -20 20 -20 25]) 三、特殊微分方程的数值求解 1、刚性微分方程的解法 振动、电路及化学反应中出现刚性现象的线性常系数方程组可表为 其中，x，f是n维向量，A是n阶矩阵。 为A的特征根。称 为刚性比，s10的方程可认为是刚性方程。 MATLAB中求解刚性方程的命令ode23s,ode15s等，其用法与ode23,ode45相同。略。