matlab将求解sin隐式解,MATLAB语言及应用-第三章课件.ppt

MATLAB语言及应用-第三章课件.ppt

第3章 MATLAB 在现代科学计算中的应用 常微分方程的MATLAB求解 一阶常微分方程的数值解法 隐式微分方程的数值解法 微分方程组的转换 刚性方程的MATLAB求解 微分代数方程求解 线性常微分方程的解析求解 一阶常微分方程的数值解法 求常微分方程数值解的函数，一般调用格式为： [t,x]=ode23(‘方程函数名’,[t0,tf], x0，选项，附加参数) [t,x]=ode45(‘方程函数名’,[t0,tf], x0，选项，附加参数) 其中t是求解的时间变量，另一个变量x为状态变量在各个时刻所组成的列向量所构成的矩阵的转置。在求解过程结束后即可以用plot(t,x)来绘制出解的结果曲线。 这里用到的方程函数名的编写格式是固定的，方程函数的引导句为 function xdot=方程函数名(t,x,flag,附加参数) 其中t为时间变量，x为状态变量，而xdot为状态变量的导数。 注意：即使微分方程是非时变的,也应该在函数输入变量列表中写上t占位。 例1 假设著名的Lorens模型的状态方程表示为 令初始条件x1(0)=x2(0)=x3(0)=10-10(1e-10),求解该微分方程。 系统动态模型： 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)]; 调用微分方程求解： t_final=100;x0=[0;0;1e-10]; [t,x]=ode45('lorenzeq', [0, t_final],x0); plot(t,x), figure; plot3(x(:,1), x(:,2), x(:,3)); 例2 考虑著名的R?ssler微分方程组 选定a=b=0.2,c=5.7,且x(0)=y(0)=z(0)=0,求解该微分方程。 隐式微分方程的数值解法 例3 假设给定的隐式微分方程如下 求解该微分方程。 原方程可以改写成矩阵形式 如果能证明矩阵A(x)为非奇异矩阵，则直接就能将该方程变换成标准的一阶微分方程组的形式 function dx=ximp(t,x) A=[sin(x(1)), cos(x(2)); -cos(x(2)), sin(x(1))]; B=[1-x(1); -x(2)];dx=inv(A)*B; -------- [t,x]=ode45(@ximp,[0,10],[0;0]); plot(t,x) 微分方程组的转换 MATLAB下提供的微分方程数值解函数只能处理一阶微分方程组形式给出的微分方程。 一阶显示微分方程的变换不是唯一的。 单个高阶常微分方程： 可以转换为： 高阶常微分方程组： 可以转换为： 例4 考虑著名的Van der Pol方程 已知y(0)=-0.2, , 求解该微分方程。 function y=vdp(t,x, mu) y=[x(2); -mu*(x(1)^2-1)*x(2)-x(1)]; global mu ----------------------- global mu x0=[-0.2;-0.7];t_final=20; mu=1; [t1,y1]=ode45(@vdp,[0,t_final],x0,mu); mu=2; [t2,y2]=ode45(@vdp,[0,t_final],x0,mu); plot(t1,y1,t2,y2, ':') figure; plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2), ':') 例5 假设系统模型以二元方程组形式给出 求解该微分方程。 刚性方程的MATLAB求解 其中一些解变化缓慢，另一些变化快，且相差较悬殊，这类方程常常称为刚性方程，又称为Stiff 方程。 刚性问题一般不适合由ode45( )这类函数求解，而应该采用ode15s( ), 调用格式和ode45( )完全一致。 例6 重新考虑著名的Van der Pol方程 已知y(0)=-0.2, , 求解该微分方程。 微分代数方程求解 一般形式 M为奇异矩阵，由求解选项中的Mass来表示该矩阵。 例如: M=[1 0 0;0 1 0;0 0 0];%质量矩阵 Options=odeset(‘mass’,m); 例6 考虑微分代数方程 并已知初始条件为x1(0)=0.8, x2(0)=x3(0)=0.1。 M=[1 0 0;0 1 0;0 0 0]; options=odeset('mass',M);