积分的数值求解函数

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,积分的数值求解函数,积分:,(求面积),trapz(x,y)给定数据点x和y，计算y=f(x)下的梯形面积积分。精度更高的求积分函数 quad,QUADL;,近似求原函数并做图mm2302,二维积分:Dblquad,三维积分:Triplequad,微分:具有其内在的敏感性,故对于被积分的数据使通过实验获得时,防止使用.应该相对数据做(),然后对()进行微分.Polyfit Polyder,diff(x)数组元素间的差值(diff(y)./diff(x)近似求微分)，结果分析mm2305,mm2306 结论,二维微分:梯度gradient 主要用作图形数据可视化(mm2307)quiver,曲率del2(mm2308),微分的数值求解函数,微分方程(ODE)1,1.IVP格式,2.IVP微分方程解法程序,t,y=ode45(,fname ,tspan,yo)用4阶/5阶龙格-库塔算法解微分方程组,other 6 ODE solvers:ODE23,ODE113,ODE15S,ODE23S,ODE23T,ODE23TB,设置选项:OPTIONS=ODESET(name1,value1,),ODEGET；,微分方程(ODE)2,output functions:,ODEPLOT,ODEPHAS2,ODEPHAS3,ODEPRINT；,evaluating solution:DEVAL；ODE,examples:RIGIDODE,BALLODE,ORBITODE,边界值问题(BVP)和偏微分方程问题,根本用法连续时间系统仿真1,1、微分方程变形：,ydot=ldefile(t,y),m文件的编写,function ydot=vdpol(t,y,mu),%VDPOL van der Pol equation.,%Ydot=VDPOL(t,Y),%Ydot(1)=Y(2),%Ydot(2)=mu*(1-Y(1)2)*Y(2)-Y(1),%mu=2,%ydot=y(2,:);mu*(1-y(1,:)2)*y(2,:)-y(1,:);,ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);,2、求解命令的编写,2401.doc,根本用法连续时间系统仿真2,例:微分方程的数值解法是在科学与工程计算中经常遇到的问题。著名的 Lorenz 模型的状态方程表示为：,假设令 ，,且假设初值为 ，,e 为一个小常数，,洛伦兹.doc,离散系统仿真,Filter,Con,SIMULINK动态仿真,SIMULINK动态仿真是MATLAB软件包之一，用于可视化的动态系统仿真，适用于连续、离散系统，也适用于线性和非线性系统。,1启动SIMULINK，进入SIMULINK窗口；,2在SIMULINK窗口下，借助模块库，创立系统 框图模型并调整模块参数。,3两种方式 SIMULINK模型窗口、matlab command 窗口sim下，设置仿真参数后，启动仿真。,4输出仿真结果：利用scope模块，利用变量返 回值、利用to workspace模块。,SIMULINK动态仿真(2),模块平衡点trim)和线性化(linmod、linmod2、dlinmod等,子系统创立及封装,S-函数,非线性系统优化设计,作业2:系统时间响应和动态仿真,应用4阶龙格-库塔法和使用simulink软件包仿真程序，求如下系统的响应。,作业2 分析求解(1),运用四阶龙格-塔库公式编程实现,公式,：,x(n+1)=x(n)+h/6(k1+2*k2+2*k3+k4),k1=fx(n),t(n),k2=fx(n)+h/2*k1,t(n)+h/2,k3=fx(n)+h/2*k2,t(n)+h/2,k4=fx(n)+h*k3,t(n)+h,编程实现,1.,定义函数dydt,function ydot=dydt(t,y),%VDPOL van der Pol equation.,%Ydot=VDPOL(t,Y),%Ydot(1)=Y(2),%Ydot(2)=Y(3),%Ydot(3)=Y(4),%Ydot(4)=Y(5),%Ydot(5)=Y(6),%Ydot(6)=,-a0*Y(1)-a1*Y(2)-a2*Y(3)-a3*Y(4)-a4*Y(5)a5*Y(6)+1,%ydot=Y(2);Y(3);Y(4);Y(5);Y(6);,-a0*Y(1)-a1*Y(2)-a2*Y(3)-a3*Y(4)-a4*Y(5)a5*Y(6)+1,编程实现,ydot=y(2);y(3);y(4);y(5);y(6);,-186.5*y(1)-840*y(2)-904.4*y(3)-237.3*y(4)-76.1*y(5)-8.8*y(6)+1,2.在workspace键入,tspan=0:0.1:5;,y0=0;0;0;0;0;0;,t,y=ode45(dydt,tspan,y0);,z=-0.2*y(:,1)+1187.6*y(:,2)+3699.6*y(:,3)+327*y(:,4)+65*y(:,5);,plot(z),或者遍写M文件 参见,work2.m,作业2 分析求解(2),用simulink系统仿真求解,1常微分方程模型时域分析,2传递函数模型,3状态方程模型,方法一：常微分方程模型,分析：,系统的输入变量为u(t),输出变量为y(t)，系统微分方程是：,D6y+8.8*D5y+76.1*D4y+237.3*D3y+904.4*D2y+840*Dy+186.5*y=65*D4u+327*D3u+3699.6*D2u+1187.6*Du-0.2*u,实现,：,运用微分模块、加法器、比例器构造系统求解,参见,work21.mdl,方法二：传递函数模型,分析：将方程两边分别取拉普拉斯变换Laplace),求传递函数G(s)。,y(t)Y(s),u(t)U(s),(a6s6+a5s5+a4s4+a3s3+a2s2+a1s+a0)Y(s),=b4s4+b3s3+b2s2+b1s+b0)U(s),G(s)=Y(s)/U(s),=b4s4+b3s3+b2s2+b1s+b0)/(a6s6+a5s5+a4s4+a3s3+a2s2+a1s+a0),实现：用传递函数模块建立系统,参见work22.mdl,方法三：状态方程模型,思路,：,复杂问题简单化。用n个一阶的微分方程去替换一个n阶的微分方程。,形式,：,X(t)=AX(t)+BU(t),Y(t)=CX(t)+DU(t),其中,X(t)=,x(1),x(2),x(n),X(t)=,x(1),x(2),x(n),=,x(2),x(3),x(n),x(n),x(n)用包含 x(1),x(2),x(n),u(t)的式子表示,状态方程模型(2),方法:重在求A、B、C、D矩阵,由常微分方程写出状态方程组,x(1)=x(2),x(2)=x(3),x(3)=x(4),x(4)=x(5),x(5)=x(6),x(6)=-a0 x(1)-a1x(2)-a2x(3)-a3x(4)-a4x(5)a5x(6)+u,y(t)=b0 x(1)+b1x(2)+b3x(3)+b4x(4),那么可以得到各矩阵,状态方程模型(3),A=,0 1 0 0 0 0,0 0 1 0 0 0,0 0 0 1 0 0,0 0 0 0 1 0,0 0 0 0 0 1,-a,0,-a,1,-a,2,-a,3,-a,4,-a,5,B=0 0 0 0 0 1,C=b,0,b,1,b,2,b,3,b,4,0,D=0,运用,simulink,的状态空间模块,输入数据求解,参见,work23.mdl,