第七章微积分的数值计算

第七章微积分的数值计算第一页，共三十三页，编辑于2023年，星期四

7.1数值微分实际问题常需计算函数的导数或积分值。但很多情况下，函数关系难以准确表示；即使能使用解析式准确表示，表示式却很复杂，不能用于实际计算。本章介绍数值计算导数或积分的实用方法。第二页，共三十三页，编辑于2023年，星期四7.1.1差分和差商

根据导数的定义其中，∆x和∆y分别称为自变量x和因变量y的增量，也称之为差分。可以用差分的商作微商（导数）的近似。数值微分就是用函数值的线性组合近似函数在某点的导数值。自变量x的步长一般取定值。首先在xi处对函数进行泰勒展开，第三页，共三十三页，编辑于2023年，星期四根据不同的组合方式可以得到精度不同的差分公式。以函数的一阶导数为例：微分公式表达式截断误差两点前向O(∆x)两点后向O(∆x)三点中心O(∆x2)三点前向O(∆x2)三点后向O(∆x2)五点中心O(∆x4)………………第四页，共三十三页，编辑于2023年，星期四精度为O(∆X2)的高阶中心差分算法精度为O(∆X4)的高阶中心差分算法第五页，共三十三页，编辑于2023年，星期四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,2……n-2。DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。第六页，共三十三页，编辑于2023年，星期四>>A=pascal(4)A=1111123413610141020>>B=diff(A)B=0123013601410>>C=diff(A,2)C=00130014>>D=diff(A,1,1)D=0123013601410>>E=diff(A,1,2)E=0001112343610第七页，共三十三页，编辑于2023年，星期四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第八页，共三十三页，编辑于2023年，星期四>>X=[13524];>>Y=gradient(X)Y=2.00002.0000-0.5000-0.50002.0000>>Y=gradient(X,2)Y=1.00001.0000-0.2500-0.25001.0000即两边用前向和后向差分，中间用中心差分>>A=pascal(4)A=1111123413610141020>>[Fx,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第九页，共三十三页，编辑于2023年，星期四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计算拉普拉斯算子。第十页，共三十三页，编辑于2023年，星期四>>4*del2(U)ans=444444444444>>[x,y]=meshgrid(1:4,1:3);>>U=x.*x+y.*yU=25101758132010131825第十一页，共三十三页，编辑于2023年，星期四7.2数值积分7.2.1数值积分基本原理

我们知道，定积分是求和式的极限，即。它的几何意义是曲边梯形的面积。从定义可知，定积分的基本分析方法是四步，即分割、近似、求和、取极限。分割就是把总量（整块曲边梯形面积）分成若干分量（小曲边梯形面积）；近似就是在每个分量中用容易计算的量去代表；求和就是把分量加起来得到总近似值；最后取极限就得到积分精确值。

第十二页，共三十三页，编辑于2023年，星期四把区间［a，b］分割成n等分，像这样取定步长算积分的方法，称为定步长积分法法。常见的低阶求积分公式复化矩形公式复化的梯形公式复化的辛普森（Simpson）公式

第十三页，共三十三页，编辑于2023年，星期四辛普森公式的几何意义

第十四页，共三十三页，编辑于2023年，星期四7.2.2变步长积分法计算积分，可以采取逐步缩小步长h的办法。即先任取步长h进行计算，然后取较小步长h’

进行计算，如果两次计算结果相差较大，则取更小步长进行计算，如此下去，直到相邻两次计算结果相差不大为止，取最小步长算出的结果作为积分值。这种方法称为变步长积分法。

利用两种步长计算积分时，通常取h’=h/2。而每次改变步长后，只需计算新增节点处的函数值，将它们的和乘新步长。

第十五页，共三十三页，编辑于2023年，星期四MATLAB常用的数值积分函数参数名功能描述quad一元函数的数值积分，采用变步长辛普森Simpson法(低阶)quadl一元函数的数值积分，采用变步长洛巴托Lobatto法(高阶)qualv一元函数的矢量数值积分dbquad二重积分（默认值为Simpson法）triplequad三重积分（默认值为Simpson法）第十六页，共三十三页，编辑于2023年，星期四

7.2.3一元函数的数值积分举例

求定积分

1.建立匿名函数

f=@(x)x./(x.^2+1);2.用sum函数实现复化矩形法求积

x=[0:0.001:1];%步长提高到0.001>>y=f(x);>>sum(y)*0.001ans=0.3468>>x=[0:0.01:1];%步长0.01

>>y=f(x);

>>sum(y)*0.01

ans=0.3491

第十七页，共三十三页，编辑于2023年，星期四3.用trapz函数实现复化梯形法求积x=[0:0.001:1];y=f(x);y=trapz(y)*0.001y=0.3466x=[0:0.01:1];y=f(x);y=trapz(y)*0.01y=0.3466第十八页，共三十三页，编辑于2023年，星期四4.用MATLAB函数求定积分用quad实现变步长辛普森法求定积分。调用格式为quad(fun,a,b,tol)其中fun为积分函数，[a,b]为积分区间，tol为积分的误差阈值，默认值为1e-6quad(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第十九页，共三十三页，编辑于2023年，星期四7.2.4矢量积分相当于多个一元定积分。例如求y=@(x,n)1./(sqrt(2*pi).*(1:n)).*exp(-x.^2./(2*(1:n).^2));%归一化高斯函数quadv(@(x)y(x,5),-1,1)ans=0.68270.38290.26110.19740.1585矢量积分的结果是一个向量，每一个元素为一个一元函数定积分的值。第二十页，共三十三页，编辑于2023年，星期四7.2.5二重定积分的数值求解使用MATLAB提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为：dblquad(f,a,b,c,d,tol,method)该函数求f(x,y)在[a,b]×[c,d]区域上的二重定积分。参数tol，可以采用method=@quadl的方法使函数用高阶的洛巴托法计算定积分。第二十一页，共三十三页，编辑于2023年，星期四例如计算f=@(x,y)exp(-x.^2-y.^2)/pi;%归一化高斯函数dblquad(f,-1,1,-1,1,1e-6,@quadl)ans=0.7101三重积分函数triplequad用法与二重积分类似第二十二页，共三十三页，编辑于2023年，星期四7.3常微分方程的数值解法7.3.1常微分方程初值问题的数值解法一般形式为所谓的数值解法就是求解y(x)在区间的近似值yn的方法，yn(n=1,2,……N)称为常微分方程的数值解。自变量x的步长一般为定值h。第二十三页，共三十三页，编辑于2023年，星期四7.3.2常见的数值方法向前欧拉公式向后欧拉公式梯形公式改进的欧拉公式第二十四页，共三十三页，编辑于2023年，星期四7.3.3龙格－库塔法简介基本思想就是利用在某些点处的值的线性组合构造公式，使其按泰勒展开后与初值问题的解的泰勒展开式比较，有尽可能多的相同项，从而保证算式有较高的精度。常用的四阶经典的龙格-库塔公式第二十五页，共三十三页，编辑于2023年，星期四7.3.4龙格－库塔法的实现求解器ODE类型特点精度说明ode45非刚性一步算法（只需前一步的结果），4,5阶Runge-Kutta方法。中大部分场合的首选算法ode23非刚性一步算法，2,3阶Runge-Kutta方法。低使用于精度较低的情形ode113非刚性多步法（需要前几布的结果），Adams算法。低~高计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法，Gear’s反向数值积分。低~中若ode45失效时，可尝试使用ode23s刚性一步法，2阶Rosebrock算法精度。低当精度较低时，计算时间比ode15s短第二十六页，共三十三页，编辑于2023年，星期四对于一个常微分方程组，如果其解相差十分悬殊，就称之为刚性方程组。对于刚性方程组，为了保持解法的稳定，步长选取十分困难，有些解法不再适用。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为初始条件。第二十七