第七章+MATLAB微积分数值计算.ppt

7.1数值微分7.2数值积分7.3常微分方程的数值解法,7.1数值微分实际问题常需计算函数的导数或积分值。但很多情况下，函数关系难以准确表示；即使能使用解析式准确表示，表示式却很复杂，不能用于实际计算。本章介绍数值计算导数或积分的实用方法。,7.1.1差分和差商根据导数的定义其中，x和y分别称为自变量x和因变量y的增量，也称之为差分。可以用差分的商作微商（导数）的近似。数值微分就是用函数值的线性组合近似函数在某点的导数值。自变量x的步长一般取定值。首先在xi处对函数进行泰勒展开，,根据不同的组合方式可以得到精度不同的差分公式。以函数的一阶导数为例：,精度为O(X2)的高阶中心差分算法,精度为O(X4)的高阶中心差分算法,7.1.2数值微分的实现在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff和梯度函数gradient。diff调用格式为：Dy=diff(Y)：计算向量Y的向前差分，并把结果赋值给向量DyDy(i)=Y(i+1)-Y(i)，i=1,2,n-1。注意向量Dy元素个数比Y少Dy=diff(Y,n)：计算向量Y的n阶向前差分。例如，diff(Y,2)=diff(diff(Y)=DX(i+1)-DX(i)=Y(i+2)-2Y(i+1)+Y(i)，i=1,2n-2。DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。差分除以设定的步长即为数值微分Dy/Dx,A=pascal(4)A=1111123413610141020B=diff(A)B=0123013601410,C=diff(A,2)C=00130014D=diff(A,1,1)D=0123013601410E=diff(A,1,2)E=0001112343610,例如：求函数f(x)=1+x2的在x110微分曲线。,见文件weifen.m,7.1.3近似梯度函数gradient的调用格式为,Fx,Fy=gradient(F,hx,hy)求矩阵F，求其x(行)方向的数值梯度Fx和y(列)方向的数值梯度Fy，x方向步长全为hx，y方向步长全为hy。Fx相当于偏导数F/x,Fy相当于偏导数F/y。Fx,Fy=gradient(F,h)求矩阵F，求其x(行)方向的数值梯度Fx和y(列)方向的数值梯度Fy，各个方向步长全为h。Fx=gradient(F)如果F是向量，直接求其数值梯度；如果F是矩阵，求其x(行)方向的数值梯度，步长为1。Fx,Fy=gradient(F)求矩阵F，求其x(行)方向的数值梯度Fx和y(列)方向的数值梯度Fy，步长为1,X=13524;Y=gradient(X)Y=2.00002.0000-0.5000-0.50002.0000Y=gradient(X,2)Y=1.00001.0000-0.2500-0.25001.0000即两边用前向和后向差分，中间用中心差分,A=pascal(4)A=1111123413610141020Fx,Fy=gradient(A)Fx=00001.00001.00001.00001.00002.00002.50003.50004.00003.00004.50008.000010.0000Fy=01.00002.00003.000001.00002.50004.500001.00003.50008.000001.00004.000010.0000,7.1.4拉普拉斯算子4*del2,由于内部算法的原因：U=4*del2(V,h)，对1维向量V以步长h求拉普拉斯算子时，返回一相同维数的向量U，且默认的步长为1。U=4*del2(V,h1,h2)，对矩阵V，横向(x方向)以步长h1，纵向(y方向)以步长h2计算拉普拉斯算子。,4*del2(U)ans=444444444444,x,y=meshgrid(1:4,1:3);U=x.*x+y.*yU=25101758132010131825,7.2数值积分7.2.1数值积分基本原理我们知道，定积分是求和式的极限，即。它的几何意义是曲边梯形的面积。从定义可知，定积分的基本分析方法是四步，即分割、近似、求和、取极限。分割就是把总量（整块曲边梯形面积）分成若干分量（小曲边梯形面积）；近似就是在每个分量中用容易计算的量去代表；求和就是把分量加起来得到总近似值；最后取极限就得到积分精确值。,把区间a，b分割成n等分，像这样取定步长算积分的方法，称为定步长积分法法。常见的低阶求积分公式复化矩形公式复化的梯形公式复化的辛普森（Simpson）公式,辛普森公式的几何意义,7.2.2变步长积分法计算积分，可以采取逐步缩小步长h的办法。即先任取步长h进行计算，然后取较小步长h进行计算，如果两次计算结果相差较大，则取更小步长进行计算，如此下去，直到相邻两次计算结果相差不大为止，取最小步长算出的结果作为积分值。这种方法称为变步长积分法。利用两种步长计算积分时，通常取h=h/2。而每次改变步长后，只需计算新增节点处的函数值，将它们的和乘新步长。,MATLAB常用的数值积分函数,7.2.3一元函数的数值积分举例求定积分,1.建立匿名函数f=(x)x./(x.2+1);2.用sum函数实现复化矩形法求积,x=0:0.001:1;%步长提高到0.001y=f(x);n=length(y);sum(y(1:n-1)*0.001ans=0.3463,x=0:0.01:1;%步长0.01y=f(x);n=length(y);sum(y(1:n-1)*0.01ans=0.3441,3.用trapz函数实现复化梯形法求积,x=0:0.001:1;y=f(x);y=trapz(y)*0.001y=0.3466,x=0:0.01:1;y=f(x);y=trapz(y)*0.01y=0.3466,4.用MATLAB函数求定积分,用quad实现变步长辛普森法求定积分。调用格式为quad(fun,a,b,tol)其中fun为积分函数，a,b为积分区间，tol为积分的误差阈值，默认值为1e-6,quad(f,0,1)ans=0.3466,用quadl实现变步长洛巴托求定积分。调用格式为quadl(fun,a,b,tol)其中fun为积分函数，a,b为积分区间，tol为积分的误差阈值，默认值为1e-6,quadl(f,0,1)ans=0.3466,7.2.4矢量积分,相当于多个一元定积分。例如求y=(x,n)1./(sqrt(2*pi).*exp(-x.2./(2*(1:n).2);%归一化高斯函数quadv(x)y(x,5),-1,1)ans=0.68270.76580.78340.78970.7926矢量积分的结果是一个向量，每一个元素为一个一元函数定积分的值。,7.2.5二重定积分的数值求解使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：dblquad(f,a,b,c,d,tol,method)该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，可以采用method=quadl的方法使函数用高阶的洛巴托法计算定积分。,例如计算f=(x,y)exp(-x.2-y.2)/pi;%归一化高斯函数dblquad(f,-1,1,-1,1,1e-6,quadl)ans=0.7101三重积分函数triplequad用法与二重积分类似,7.3常微分方程的数值解法,7.3.1常微分方程初值问题的数值解法一般形式为所谓的数值解法就是求解y(x)在区间的近似值y的方法，y(n=1,2,N)称为常微分方程的数值解。自变量x的步长一般为定值h。,7.3.2常见的数值方法向前欧拉公式向后欧拉公式梯形公式改进的欧拉公式,7.3.3龙格库塔法简介基本思想就是利用在某些点处的值的线性组合构造公式，使其按泰勒展开后与初值问题的解的泰勒展开式比较，有尽可能多的相同项，从而保证算式有较高的精度。常用的四阶经典的龙格-库塔公式,7.3.4龙格库塔法的实现,对于一个常微分方程组，如果其解相差十分悬殊，就称之为刚性方程组。对于刚性方程组，为了保持解法的稳定，步长选取十分困难，有些解法不再适用。ode函数调用格式：t,y=odeij(odefun,tspan,y0)t,y=odeij(odefun,tspan,y0,options)t,y=odeij(odefun,tspan,y0,options,p1,p2,.)t,y,te,ye,ie=odeij(odefun,tspan,y0,options,p1,p2,.)odefun为显式常微分方程中的f(xn,yn)，tspan为求解区间，y0为初始条件。,7.3.5求解多变量一阶常微分方程组。例：求解下面的混沌理论的洛仑兹方程。x(0)=0;y(0)=0;z(0)=0;,编写lorfun.m如下：functionydot=lorfun(t,y)ydot=-8/3*y(1)+y(2)*y(3);-10*y(2)+10*y(3);-y(2)*y(1)+35*y(2)-y(3);建立lorfun.m文件t,y=ode23(lorfun,0,20,0,0,eps);plot3(y(:,1),y(:,2),y(:,3)gridon,7.3.6求解高阶常微分方程可以先把其转化为一阶常微分方程组,例令则,首先编写vdp.m如下：functionfy=vdp(x,y