数值积分与数值微分.docx

1、数值积分与数值微分8 数值积分与数值微分8.1 例题解答例 8.1 给定积分,分别用梯形公式、公式、公式作近似计算.解:先输入主要初始参数a=0.5;b=1;f=inline(x(1/2);%梯形公式I1=(b-a)/2*(feval(f,a)+feval(f,b)I1 = 0.426776695296637%simpson公式I2=(b-a)/6*(feval(f,a)+4*feval(f,(a+b)/2)+feval(f,b)I2 = 0.430934033027025%Cotes公式(n=4)tc=0;C0=7 32 12 32 7;for i=0:4 tc=tc+C0(i+1)*fev

2、al(f,a+i*(b-a)/4);endI3=(b-a)/90*tc I3 = 0.430964070495876%准确值I=int(char(f),a,b)vpa(I) I =-1/6*2(1/2)+2/3ans =0.43096440627115082519971854596505例 8.2 对积分,为使其精度达到.若用复化梯形公式,应将0,1多少等分?若用复化公式,应将0,1多少等分?解:直接按余项计算即可.复化梯形公式的余项为:复化公式余项为:对于,在课本中我们已证得以下不等式成立:直接利用上述不等式关系解答本题.先输入误差精度:Eps=1E-4 Eps = 1.0000000000

3、00000e-004(1) 复化梯形公式h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1) %先求出步长h1 = 0.060000000000000 N1=ceil(1/h1) %向上取整,得到等分区间数N1 = 17故可将区间17等分即可达到所要求的精度.(2) 复化公式h2=power(Eps/abs(-(1-0)/180*1/(1+4),1/4) %先求出步长h2 = 0.547722557505166 N2=ceil(1/h2) %向上取整,得到等分区间数N2 = 2 故可将区间2等分即可达到所要求的精度. 扩展:1) Matlab中复化梯形公式命令为I=trapz(x

4、,y),复化公式命令为quad().2) Matlab中有四个取整函数,分别为ceil(),floor(),fix(),round(),分别表示向正无穷大方向取整、向负无穷大方向取整、向靠近零方向舍入和四舍五入.例 8.3 对积分,利用变步长方法求其近似值,使其精度达到.解：利用变步长法前先建立三种变步长复化积分公式的函数. 注意在Matlab中直接用sin(0)/0得不到1，因此解此题时我们改用求极限的方法得到函数值,此函数名为limit().先建立三种复化公式的函数文件,它们分别为复化梯形公式trap.m、复化公式为simpson.m、公式为cotes.m,三个函数的源程序如下:(1) 复

5、化梯形公式trap.mfunction T=trap(f,a,b,n)%trap.m%复化梯形公式求积分值%f为积分函数%a,b为积分区间%n是等分区间份数h=(b-a)/n;%步长T=0;for k=1:(n-1) x0=a+h*k; T=T+limit(f,x0);endT=h*(limit(f,a)+limit(f,b)/2+h*T;T=double(T);(2) 复化公式simpson.m:function S=simpson(f,a,b,n)%simpson.m%Simpson公式求积分值%f为积分函数%a,b为积分区间%n是等分区间份数h=(b-a)/(2*n);%步长s1=0;s

6、2=0;for k=1:n x0=a+h*(2*k-1); s1=s1+limit(f,x0);endfor k=1:(n-1) x0=a+h*2*k; s2=s2+limit(f,x0);endS=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3;S=double(S);(3) 复化公式cotes.m:function C=cote(f,a,b,n)%cote.m%复化cotes公式求积分值%f为积分函数%a,b为积分区间%n是等分区间份数h=(b-a)/n;%步长C=0;for i=1:(n-1) x0=a+i*h; C=C+14*limit(f,x0);endf

7、or k=0:(n-1) x0=a+h*k; s=32*limit(f,x0+h*1/4)+12*limit(f,x0+h*1/2)+32*limit(f,x0+h*3/4); C=C+s;endC=C+7*(limit(f,a)+limit(f,b);C=C*h/90;C=double(C);再编写主程序调用这三个函数,主程序名为ex8_3.m,源程序如下:%ex8_3.mclc;syms x;f=sym(sin(x)/x);a=0;b=1;%积分上下限n=20;%作1,2,3,20次区间等分%复化梯形公式T=zeros(n,1);for i=1:n T(i)=trap(f,a,b,i);e

8、nd%复化Simpson公式;S=zeros(n,1);for i=1:n S(i)=simpson(f,a,b,i);end%复化Cotes公式C=zeros(n,1);for i=1:n C(i)=cote(f,a,b,i);end%准确值I=int(f,a,b);I=double(I);%画图作出直观观察x=;x=1:n;figure;plot(x,ones(1,n)*I,-);hold on;plot(x,T,r-,LineWidth,2);plot(x,S,m.-,LineWidth,1);plot(x,C,c:,LineWidth,1.5);grid on;title(三种复化公式

9、积分效果对比图);legend(准确值曲线,复化梯形公式,复化Simpson公式,复化Cotes公式);hold off;disp( 复化梯形公式 复化Simpson公式 复化Cotes公式);disp(T S C);在Matlab命令窗口中输入ex8_3,得到如下积分结果:复化梯形公式 复化Simpson公式 复化Cotes公式 0.920735492403948 0.946145882273587 0.946083004063674 0.939793284806177 0.946086933951794 0.946083069350917 0.943291429132337 0.94608

10、3831311699 0.946083070278278 0.944513521665390 0.946083310888472 0.946083070351379 0.945078780953402 0.946083168838073 0.946083070363043 0.945385730766859 0.946083117842867 0.946083070365797 0.945570776256246 0.946083095989403 0.946083070366633 0.945690863582701 0.946083085384948 0.946083070366936 0

11、.945773188549752 0.946083079742053 0.946083070367061 0.945832071866905 0.946083076517732 0.946083070367118 0.945875637119198 0.946083074567937 0.946083070367147 0.945908771073865 0.946083073333114 0.946083070367161 0.945934556488475 0.946083072520476 0.946083070367170 0.945955016056114 0.94608307196

12、8057 0.946083070367174 0.945971521552985 0.946083071581964 0.946083070367177 0.945985029934386 0.946083071305562 0.946083070367179 0.945996225242376 0.946083071103489 0.946083070367180 0.946005606943978 0.946083070952998 0.946083070367181 0.946013546623966 0.946083070839065 0.946083070367182 0.94602

13、0325355025 0.946083070751532 0.946083070367182由以上结果可看到复化梯形公式有一个上升接近准确值过程,而复化公式积分结果和复化公式积分的结果基本上和准确值的曲线重叠在一块,可见它们的精度是相当高的.例 8.4 用积分法计算,精度.解：编写积分法的函数M文件,源程序如下(romberg.m):function I,T=romberg(f,a,b,n,Eps)%Romberg积分计算%f为积分函数%a,b为积分区间%n+1是T数表的列数目%Eps为迭代精度%返回值中I为积分结果,T是积分表if narginEps) & (jn)| (jsyms x;%创

14、建符号变量f=sym(sin(x)/x) %符号函数f =sin(x)/x I,T=romberg(f,0,1,3,1E-6) %积分计算 I = 0.9461T = 0.9207 0 0 0 0 0.9398 0.9461 0 0 0 0.9445 0.9461 0.9461 0 0 0.9457 0.9461 0.9461 0.9461 0 0.9460 0.9461 0.9461 0.9461 0.9461 其中T为积分表,由输出结果可知计算结果为.例 8.5 利用积分法求.解:直接利用例8.4的函数即可.syms x;f=sym(4/(1+x2) f =4/(1+x2) I,T=rom

15、berg(f,0,1,5,1E-6) I = 3.1416T = 3.0000 0 0 0 0 0 3.1000 3.1333 0 0 0 0 3.1312 3.1416 3.1421 0 0 0 3.1390 3.1416 3.1416 3.1416 0 0 3.1409 3.1416 3.1416 3.1416 3.1416 0 3.1414 3.1416 3.1416 3.1416 3.1416 3.1416 故结果为.例 8.6 对积分构造其型求积公式.解:取,得到方程组:这是一个抽象方程组,可以利用Matlab的符号法来解之,该函数名为solve():%直接输入解之:x=solve(

16、A0+A1=2,A0*x0+A1*x1=0,A0*x02+A1*x12=2/3,A0*x03+A1*x13=0,A0,A1,x0,x1) x = A0: 2x1 sym A1: 2x1 sym x0: 2x1 sym x1: 2x1 sym%显示结果x.A0,x.A1,x.x0,x.x1 ans = 1 1ans = 1 1ans = 1/3*3(1/2) -1/3*3(1/2)ans = -1/3*3(1/2) 1/3*3(1/2)因此得到两组解为:或: 求积公式为例 8.7 分别利用公式及公式计算积分:解:(1)先计算准确值: I=int(x+1.5)(1/2),-1,1) I =2.39

17、95291230781336018654631663256(2)两点公式:由课本表中查得两点公式的求积系数为1,求积节点为f=inline(x+1.5)(1/2)f = Inline function: f(x) = (x+1.5)(1/2) I1=1*feval(f,-0.57735)+1*feval(f,0.57735) I1 = 2.401848214499055(3)两个节点梯形公式:f=inline(x+1.5)(1/2)a=-1;%求积区间b=1;I2=(b-a)/2*(feval(f,a)+feval(f,b) %两点梯形公式f = Inline function: f(x) =

18、 (x+1.5)(1/2)I2 = 2.288245611270737(4)三点公式:查表有,:f=inline(x+1.5)(1/2)x=0.7745966692,- 0.7745966692,0%节点系数A=0.555555556, 0.555555556,0.888888889%求积系数I3=0;for i=1:length(x) I3=I3+feval(f,x(i)*A(i);%三点Gauss-Legendre公式endI3 %显示结果f = Inline function: f(x) = (x+1.5)(1/2)x = 0.774596669200000 -0.77459666920

19、0000 0A = 0.555555556000000 0.555555556000000 0.888888889000000I3 = 2.399708072133707 (5)三个节点公式:f=inline(x+1.5)(1/2)a=-1;%积分区间b=1;x=a,(a+b)/2,b%积分节点A=1 4 1%积分节点系数I4=0;for i=1:length(x) I4=I4+feval(f,x(i)*A(i);%三点Simpson公式endI4=I4*(b-a)/6 %显示计算结果f = Inline function: f(x) = (x+1.5)(1/2)x = -1 0 1A = 1

20、 4 1I4 = 2.395741698945698例 8.8 利用求积公式求积分.(准确值为1.)解: 先作变换,将积分区间变换到上,令,则有:于是就可以利用求积公式来解题.(1) 两点求积公式: f=inline(sin(pi*(t+1)/4) f = Inline function: f(t) = sin(pi*(t+1)/4) I1=1*feval(f,-0.57735)+1*feval(f,0.57735);I1=I1*pi/4 I1 = 0.998472716275797(2)三点求积公式:f=inline(sin(pi*(t+1)/4) f = Inline function:

21、f(t) = sin(pi*(t+1)/4) x=0.7745966692,- 0.7745966692,0%节点系数 x = 0.774596669200000 -0.774596669200000 0 A=0.555555556, 0.555555556,0.888888889%求积系数 A = 0.555555556000000 0.555555556000000 0.888888889000000 I2=0;for i=1:length(x)I2=I2+feval(f,x(i)*pi/4*A(i);endI2 I2 = 1.000008122033779例 8.9 计算积分,精度要求.

22、解:先利用余项来估计误差精度,再进行计算即可.这里利用求积公式余项作误差估计.余项为:对右边进行放大得到:%先搜索满足上述不等式的项数ni=1;E= pi*exp(1)/(power(2,2*i+1)*factorial(2*i+2);while E1E-6 i=i+1;E= pi*exp(1)/(power(2,2*i+1)*factorial(2*i+2);endi,E %显示结果i = 4E = 4.596331680902589e-009可见只要计算前4项即可.f=inline(exp(x) f = Inline function: f(x) = exp(x) n=i;I=0;for j=0:n x=cos(pi*(2*j+1)/(2*n+2); I=I+feval(f,x);endI=I*pi/(n+1) I = 3.977463258776695