《第五章 数值积分-席.ppt》由会员分享，可在线阅读，更多相关《第五章 数值积分-席.ppt（57页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、对定积分进行对定积分进行MATLAB符号计算可以使用命令：符号计算可以使用命令：f=int(s,a,b)输入量：输入量：s可以是被积函数，也可以是几个被积可以是被积函数，也可以是几个被积函数组成的矩阵；积分上下限函数组成的矩阵；积分上下限a，b可以是数也可以是可以是数也可以是符号变量；积分变量默认为符号变量；积分变量默认为x，若积分变量为，若积分变量为t，采用，采用f=int(s,t,a,b)第五章第五章 数值积分数值积分例例:计算计算解：解：syms t f=5/(t-1)*(t-2)*(t-3);F=int(f,t,4,5),y=double(F)例：例：解：解：syms x a b c
2、d nf=(a*x+b)/(c*x+d)(1/n);F=int(f,0,1),y=simple(F)（1）积不出：如）积不出：如 等；等；在求定积分在求定积分 时，时，（2）难积：当）难积：当 比较复杂时，不好使用比较复杂时，不好使用NL公式；公式；（3）表格函数）表格函数问题问题：一、梯形法一、梯形法对积分对积分用线性插值函数用线性插值函数 近似代替被积函数近似代替被积函数 ，其中：其中：则则上式称为上式称为梯形积分公式。梯形积分公式。也可写为：也可写为：其中其中复合梯形积分法复合梯形积分法：误差误差 若被积函数在等距的若被积函数在等距的 n1 个点上取值，可个点上取值，可在每个区间段逐个利
3、用梯形公式，整理后得：在每个区间段逐个利用梯形公式，整理后得：Matlab里提供了复合梯形公式的计算程序里提供了复合梯形公式的计算程序trapz.m，其调用格式为：，其调用格式为：Z=trapz(y）y为向量，则为向量，则Z是是y的积分，若的积分，若y为矩阵，则行向量为矩阵，则行向量Z是是y的每一列的积分；的每一列的积分；Z=trapz(x,y)输入输入x、y为同长度的数组，计算为同长度的数组，计算y对对x积分。积分。例例：X=0:pi/100:pi;Y=sin(X);Z=trapz(X,Y)Z0=pi/100*trapz(Y)Z=1.9998Z0=1.9998Y=1 2 3 4;trapz(
4、Y)ans=7.5000根据上面的复合梯形积分公式：trapz(Y)=1+2(2+3)+4/2基于梯形法的积分函数文件有基于梯形法的积分函数文件有trapez_v，trapez_n和和trapez_g，分别如下：，分别如下：function I=trapez_n(f_name,a,b,n)h=(b-a)/n;x=a+(0:n)*h;f=feval(f_name,x);%求求f_name(x),f_name被积函数被积函数%文件名文件名:feval(fun,x)或或feval(fun,x)I=trapez_v(f,h)function I=trapez_v(f,h)I=h*(sum(f)-(f(
5、1)+f(length(f)/2);%复合梯形公式复合梯形公式function I=trapez_g(f_name,a,b,n)h=(b-a)/n;x=a+(0:n)*h;f=feval(f_name,x);I=h/2*(f(1)+f(n+1);%n=1,积分区间没有进行分割积分区间没有进行分割 if n1,I=I+h*sum(f(2:n);end I;h2=(b-a)/100;xc=a+(0:100)*h2;fc=feval(f_name,xc);plot(xc,fc,r);hold on title(Trapezoidal Rule);xlabel(x);ylabel(y);plot(x,
6、f,ko);plot(x,zeros(size(x)%x轴上分点图形轴上分点图形 for i=1:n;plot(x(i),x(i),0,f(i);end例例1 已知积分精确值已知积分精确值 ，分析用梯形法计算，分析用梯形法计算 下述积分时，分割区间数下述积分时，分割区间数 n 对误差的影响。对误差的影响。clear;Iexact=4.006994;a=0;b=2;n=1;fprintf(n n I Errorn);for k=1:6 n=2*n;h=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=sqrt(1+exp(x);I=trapez_v(f,h);fprintf(%3.0f
7、%10.5f%10.5fn,n,I,Iexact-I);end运行结果：运行结果：n I Error 2 4.08358 -0.07659 4 4.02619 -0.01919 8 4.01180 -0.00480 16 4.00819 -0.00120 32 4.00729 -0.00030 64 4.00707 -0.00008结结果果表表明明：复复合合梯梯形形公公式式的的误误差差随随着着 n 的的加加倍倍减减少少为为原来的原来的1/4。调用trapez_g函数：第一步建立被积函数文件：函数：第一步建立被积函数文件：function y=f(x)y=sqrt(1+exp(x)再调用：再调用
8、：I=trapez_g(f,0,2,16)上例表明：复合梯形公式的误差随着上例表明：复合梯形公式的误差随着 n 的加倍的加倍减少为原来的减少为原来的1/4，可以利用这种趋势来降低主要误，可以利用这种趋势来降低主要误差部分。记差部分。记 为分割区间数为为分割区间数为 n 的积分值，的积分值，为区为区间加倍后的积分值，则间加倍后的积分值，则 约为约为 误差的三倍，误差的三倍，减掉这一误差，结果将会更加精确。减掉这一误差，结果将会更加精确。这种方法称为这种方法称为龙贝格积分龙贝格积分，记为：，记为：由此而得的结果相当于由此而得的结果相当于 2n 个区间上的个区间上的1/3辛普森法辛普森法。二、辛普森
9、法二、辛普森法1、1/3辛普森法辛普森法对积分对积分 其中其中可得：可得：若用若用二次多项式二次多项式代替代替 若加误差项若加误差项称为称为1/3辛普森法辛普森法。若将积分区间分为若将积分区间分为 n（偶数）个小区间，在每个小（偶数）个小区间，在每个小区间上逐个用区间上逐个用1/3辛普森法，整理得：辛普森法，整理得：其中其中复合复合1/3辛普森积分公式辛普森积分公式例例2 用复合用复合1/3辛普森求积公式计算辛普森求积公式计算clearIexact=4.006994;a=0;b=2;n=1;fprintf(n I Errorn);for k=1:4 n=2*n;h=(b-a)/n;i=1:n+
10、1;x=a+(i-1)*h;f=sqrt(1+exp(x);I=(h/3)*(f(1)+4*sum(f(2:2:n)+f(n+1);if n2,I=I+(h/3)*2*sum(f(3:2:n);end fprintf(%3.0f%10.5f%10.5fn,n,I,Iexact-I);end运行结果运行结果 n I Error 2 4.00791 -0.00092 4 4.00705 -0.00006 8 4.00700 -0.00000 16 4.00699 -0.000002、3/8辛普森法辛普森法若对积分，若对积分，其中其中3/8辛普森法辛普森法类似可得类似可得复合复合3/8辛普森法辛普森
11、法三等分积分区间，构造出三等分积分区间，构造出三次插值多项式三次插值多项式代替代替 ，则，则注注：对辛普森法，可建立函数对辛普森法，可建立函数 simps_v 和和 simps_n 实现。实现。function I=Simps_n(f_name,a,b,n)h=(b-a)/n;x=a+(0:n)*h;f=feval(f_name,x);I=Simps_v(f,h)(1)复合复合3/8辛普森法仅适用于辛普森法仅适用于分割区间为分割区间为3的倍的倍 数数的情形，一般看做的情形，一般看做1/3法的补方；法的补方；(2)当分割区间为奇数时，可当分割区间为奇数时，可 将前三个或后三个区将前三个或后三个区
12、 间用间用3/8法，其余区间仍用法，其余区间仍用1/3法。两种方法结合法。两种方法结合 使用不降低积分的精确度；使用不降低积分的精确度；function I=Simps_v(f,h)n=length(f)-1;if n=1,fprintf(Data has only one interval),return;endif n=2,I=h/3*(f(1)+4*f(2)+f(3);return;endif n=3,I=3/8*h*(f(1)+3*f(2)+3*f(3)+f(4);return;endI=0;if 2*floor(n/2)=n I=3/8*h*(f(n-2)+3*f(n-1)+3*f(
13、n)+f(n+1);m=n-3else m=nendI=I+(h/3)*(f(1)+4*sum(f(2:2:m)+f(m+1);if m2,I=I+(h/3)*2*sum(f(3:2:m);end%如果奇数等分，如果奇数等分，%则后三个区间采用则后三个区间采用3/8，其余采用，其余采用1/3辛普森积分公式辛普森积分公式例例3 计算积分计算积分R=5;g=9.81;b=0.1;z1=-0.9*R;z2=R;h=(z2-z1)/20;z=z1:h:z2;f=(R2-z.2)./(b2*sqrt(2*g*(z+R);I=Simps_v(f,h)结果结果三、三、Newton-Cotes公式公式 New
14、ton-Cotes公式是由拉格朗日插值公式而推公式是由拉格朗日插值公式而推导出来的一个系列的数值积分公式。导出来的一个系列的数值积分公式。将区间将区间 分为分为 n 等份，记等份，记 ，分点为，分点为,这这 n+1个节点个节点 上的函上的函数值为数值为 ,从而区间从而区间 上的拉格朗日插值多项式为上的拉格朗日插值多项式为其中，其中，为插值基函数，与函数为插值基函数，与函数 无关。无关。其中其中以以 近似被积函数近似被积函数 得：得：设设则则对对将其被积函数转换为幂级数形式：将其被积函数转换为幂级数形式：上式左边是多项式，上式左边是多项式，所以系数所以系数 由由polyfit函数确定。函数确定。
15、s=1:n+1;y=eye(n+1);for i=1:n+1a(i,:)=polyfit(s,y(i,:),n);enda满足：满足：因此因此此式称为此式称为Newton-Cotes求积公式求积公式，称为权系数。称为权系数。（称为称为Cotes系数）系数）由上述过程计算出来的由上述过程计算出来的Cotes系数及余项如下表：系数及余项如下表：5 6余项余项Cotes系数系数 n当当 n 4时，时，NewtonCotes公式为：公式为：Cotes公式NewtonCotes积分公式可以由函数积分公式可以由函数Newt_itg实现：实现：function I=Newt_itg(f_name,a,b,n
16、)npt=n+1;en=npt:-1:1;x=1:npt;for i=1:npt power_2(i)=npten(i);power_1(i)=1en(i);endfor j=1:npt z=zeros(1,npt);z(j)=1;a1=polyfit(x,z,npt-1);w(j)=sum(a1.*(power_2-power_1)./en);endx=a:(b-a)/(npt-1):b;y=feval(f_name,x);I=sum(w.*y)*(b-a)/(npt-1);fprintf(n x y w n)for j=1:nptfprintf(%e%e%en,x(j),y(j),w(j)
17、end 1.建立被积函数文件：建立被积函数文件：function y=f(x)y=sqrt(1+exp(x)2.调用：调用：I=Newt_itg(f,0,2,5)例：四、四、GaussLegendre公式公式 采用采用不等距节点不等距节点（Legendre多项式的根）可以多项式的根）可以大幅度提高积分精度，在相同的节点数下，其精度大幅度提高积分精度，在相同的节点数下，其精度是是NewtonCotes公式的两倍。公式的两倍。Legendre多项式：多项式：function pn=Legen_pw(n)pbb=1;if n=0,pn=pbb;break;endpb=1 0;if n=1,pn=pb
18、;break;endfor i=2:n;pn=(2*i-1)*pb,0-(i-1)*0,0,pbb)/i;pbb=pb;pb=pn;end调用：p=Legen_pw(n)其幂级数形式为：其幂级数形式为：（系数可由（系数可由Legen_pw(n)确定）确定）n 阶阶Legendre多项式的多项式的 n 个根全部位于个根全部位于-1与与1之之间，间，n7时的根如下，在端点处没有根，因此向两时的根如下，在端点处没有根，因此向两端点方向区间长度逐渐减少。端点方向区间长度逐渐减少。对积分对积分用用 n 个个Legendre节点的插值多项式拟合被积函数，节点的插值多项式拟合被积函数，然后积分：然后积分：其
19、中，其中，是是 n 阶阶Legendre多项式的根，多项式的根，满足：满足：而系数而系数 由由polyfit函数确定。函数确定。当积分限不是当积分限不是-1和和1时，通过坐标变换可将这些根映时，通过坐标变换可将这些根映射到积分上、下限射到积分上、下限 a，b之间，如：之间，如：Legendre求积公式：求积公式：Legendre多项式的根可由多项式的根可由roots(Legen_pw(n)给出，给出，该求积公式由函数该求积公式由函数Gauss_q实现。实现。function I=Gauss_q(f_name,a,b,n)p=legen_pw(n);x=roots(p);x=sort(x);fo
20、r j=1:n y=zeros(1,n);y(j)=1;p=polyfit(x,y,n-1);P=poly_itg(p);%求原函数求原函数 w(j)=polyval(P,1)-polyval(P,-1);endx=0.5*(b-a)*x+a+b);y=feval(f_name,x);I=sum(w.*y)*(b-a)/2;fprintf(n x y w n)for j=1:n，fprintf(%e%e%en,x(j),y(j),w(j)；endfunction py=poly_itg(p)n=length(p);py=p.*n:-1:1.(-1),0例例4：用：用6个节点的个节点的Gauss
21、公式计算积分公式计算积分I=Gauss_q(fex,0,pi,6)function f=fex(x)f=1./(2+cos(x);结果：结果：I1.8138五、积分限无界及被积函数有奇点的数值积分法五、积分限无界及被积函数有奇点的数值积分法考虑如下积分考虑如下积分无穷区间的积分无穷区间的积分无界函数的积分无界函数的积分 一个函数在无穷区间或半无穷区间可积，当且一个函数在无穷区间或半无穷区间可积，当且仅当该函数在仅当该函数在0附近的小区域内取值大幅度变化，而附近的小区域内取值大幅度变化，而当当 x 趋于无穷（负无穷）时，函数值趋于趋于无穷（负无穷）时，函数值趋于0。先计算先计算其中其中 X 足够
22、大使得位于足够大使得位于 之外的区间对于之外的区间对于积分值的影响可以忽略不计。积分值的影响可以忽略不计。计算积分计算积分1、复合梯形求积公式的使用、复合梯形求积公式的使用-无穷区间的积分无穷区间的积分对积分对积分最有效的积分方法是复合梯形公式：最有效的积分方法是复合梯形公式：其中其中复合梯形公式使用较少的节点达到相对较高的精度。复合梯形公式使用较少的节点达到相对较高的精度。例例5 计算积分计算积分解：解：调用函数调用函数trapez_v，分别取分别取 n10和和20，程序如下：，程序如下：用用 10 和和 10替换积分限得：替换积分限得：function I=trapez_v(f,h)I=h
23、*(sum(f)-(f(1)+f(length(f)/2);%复合梯形公式复合梯形公式cleara=-10;b=10;n=5;fprintf(n h n I n);for k=1:2 n=2*n;h=(b-a)/n;i=1:n+1;x=a+(i-1)*h;f=exp(-x.2)/sqrt(pi);I=trapez_v(f,h);fprintf(%3.0f%10.5f%10.5fn,h,n,I);endh n I2 10 1.169711 20 1.00010 注：通常，注：通常，X 的取值不能预先知道，而是通过一些的取值不能预先知道，而是通过一些试验的方法取得的。试验的方法取得的。X 的最优取
24、值应使得，当的最优取值应使得，当 X 减减少少1.5倍时，倍时，I 值会受到影响，而增大值会受到影响，而增大1.5倍时，倍时，I 值值不会发生变化。不会发生变化。2、指数变换、指数变换-无界函数的积分无界函数的积分 当被积函数在积分限的一端或两端有奇点时，可当被积函数在积分限的一端或两端有奇点时，可将有限区域将有限区域 a,b 通过坐标变换转换为无穷区间的积通过坐标变换转换为无穷区间的积分，再利用复合梯形求积公式。分，再利用复合梯形求积公式。对积分对积分设映射设映射使得使得则则如设如设或或则则例例6：计算：计算令令取等步长取等步长 h0.2，分割区间，分割区间 ，节点值为：，节点值为：其中其中
25、则则下图显示了下图显示了 z 坐标下等距节点到坐标下等距节点到 x 坐标的转换坐标的转换:下图显示了被积函数的节点及其图形下图显示了被积函数的节点及其图形:下图显示了梯形公式的节点及其图形下图显示了梯形公式的节点及其图形:六、六、MATLABMATLAB中的积分命令中的积分命令 MATLAB工具箱中的积分命令有工具箱中的积分命令有quad和和quad8，前者采用循环辛普森法，后者采用循环八阶，前者采用循环辛普森法，后者采用循环八阶NewtonCotes法，调用格式为：法，调用格式为：quad(f_name,a,b)quad(f_name,a,b,tol)quad(f_name,a,b,tol,
26、trace)第一种格式中，容许误差设为默认值第一种格式中，容许误差设为默认值 0.000001，求，求积运算一直进行到容许误差满足为止。如果第三种积运算一直进行到容许误差满足为止。如果第三种格式中取格式中取trace非非0，则积分进程将显示在屏幕上。，则积分进程将显示在屏幕上。quad8的调用格式与的调用格式与quad基本相同。基本相同。注：函数注：函数quad(f_name,a,b,tol,trace)用于求被积函用于求被积函数数f(x)在在a,b上的定积分，上的定积分，tol是计算精度，调用是计算精度，调用quad函数时，先要建立一个描述被积函数函数时，先要建立一个描述被积函数f(x)的函
27、的函数文件或语句函数。当被积函数数文件或语句函数。当被积函数f_name含有一个以含有一个以上的变量时，上的变量时，quad函数的调用格式为：函数的调用格式为：quad(f_name,a,b,tol,trace,g1,g2)其中其中f_name,a,b,tol,trace等参数的含义同前。等参数的含义同前。例例 用两种不同的方法求积分。用两种不同的方法求积分。先建立一个函数文件先建立一个函数文件ex.m：function ex=ex(x)ex=exp(-x.2);%注意应用点运算注意应用点运算return或者直接用或者直接用ex1=inline(exp(-x.2);然后，在然后，在MATLAB
28、命令窗口，输入命令：命令窗口，输入命令：quad(ex,0,1,1e-6)%注意函数名应加字符引号注意函数名应加字符引号或者或者quad(ex1,0,1,1e-6)quad8(ex,0,1,1e-6)%用另一函数求积分用另一函数求积分注：注：inline函数可以直接定义函数的表达式，不需函数可以直接定义函数的表达式，不需要使用函数的要使用函数的m文件。文件。另外，另外，MATLAB还提供了计算二重积分的函数：还提供了计算二重积分的函数：dblquad(f,a,b,c,d,tol,trace)该函数求该函数求f(x,y)在在a,bc,d区域上的二重积分。参区域上的二重积分。参 数数tol，tra
29、ce的用法与函数的用法与函数quad完全相同。完全相同。命令如下：命令如下：g=inline(exp(-x.2-y.2);dblquad(g,0,1,0,1)%直接调用二重积分函数求解直接调用二重积分函数求解h=0.2;i=-30:30;zi=i*h;xi=(1+tanh(zi)/2;plot(xi,zi);hold onfor j=1:length(i)plot(xi(j),xi(j),-6,zi(j);plot(0,xi(j),zi(j),zi(j);endhold offtitle(x_i=(1+tanh(z_i)/2);xlabel(x);ylabel(z);z 坐标下等距节点到坐标下
30、等距节点到 x 坐标的转换坐标的转换:h=0.2;i=-30:30;zi=i*h;xi=(1+tanh(zi)/2;y=exp(-xi.2)./sqrt(1-xi.2);plot(xi,y,-,xi,y,ro)axis(0,1,0,7);legend(被积函数的节点及其图形)xlabel(x,Color,w);ylabel(y,Color,w);被积函数的节点及其图形被积函数的节点及其图形:h=0.2;i=-30:30;zi=i*h;xi=(1+tanh(zi)/2;y=exp(-xi.2)./sqrt(1-xi.2);g_z=y./(2*(cosh(zi).2);plot(zi,g_z,-,zi,g_z,ro)axis(-6,6,0,0.7);legend(梯形公式的节点及其图形)xlabel(z,Color,w);ylabel(g(z),Color,w);梯形公式的节点及其图形梯形公式的节点及其图形: