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)*feval(f,a+i*(b-a)/4); end
>>I3=(b-a)/90*tc I3 =
0.430964070495876
%准确值
>>I=int(char(f),a,b) >>vpa(I) I =
-1/6*2^(1/2)+2/3 ans =
0.43096440627115082519971854596505
?10.5xdx,分别用梯形公式、Simpson公式、Cote公式作近似计算.
例 8.2 对积分I?sinx?4?0xdx,为使其精度达到10.
1若用复化梯形公式,应将[0,1]多少等分?
若用复化Simpson公式,应将[0,1]多少等分? 解:
直接按余项计算即可. 复化梯形公式的余项为:
En(f)??复化Simpson公式余项为:
b?a2''hf(?),12??(a,b)
b?a?h?(4)En(f)????f(?),180?2?对于f(x)?4??(a,b),f(x)?C4[a,b]
sinx,在课本中我们已证得以下不等式成立: x1f(k)(x)?
k?1直接利用上述不等式关系解答本题. 先输入误差精度: >>Eps=1E-4
Eps =
1.000000000000000e-004
(1) 复化梯形公式
>>h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1))) %先求出步长
h1 =
0.060000000000000
>>N1=ceil(1/h1) %向上取整,得到等分区间数 N1 = 17
故可将区间17等分即可达到所要求的精度. (2) 复化Simpson公式
>>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,y),复化Simpson公式命令为quad().
2) Matlab中有四个取整函数,分别为ceil(),floor(),fix(),round(),分别表示向正无穷大方向取
整、向负无穷大方向取整、向靠近零方向舍入和四舍五入.
例 8.3 对积分I?sinx?6dx,利用变步长方法求其近似值,使其精度达到. ??10?0x1解：
利用变步长法前先建立三种变步长复化积分公式的函数. 注意在Matlab中直接用sin(0)/0得不到1，此函数名为limit().
先建立三种复化公式的函数文件,它们分别为复化梯形公式trap.m、复化Simpson公式为simpson.m、Cotes公式为cotes.m,三个函数的源程序如下: (1) 复化梯形公式trap.m
function 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); end
T=h*(limit(f,a)+limit(f,b))/2+h*T; T=double(T);
(2) 复化Simpson公式simpson.m:
function S=simpson(f,a,b,n) %simpson.m
%Simpson公式求积分值 %f为积分函数 %[a,b]为积分区间 %n是等分区间份数
h=(b-a)/(2*n);%步长 s1=0; s2=0; for k=1:n
sin(0)0matlab?NAN，因此解此题时我们改用求极限的方法得到函数值,
x0=a+h*(2*k-1); s1=s1+limit(f,x0); end
for k=1:(n-1)
x0=a+h*2*k; s2=s2+limit(f,x0); end
S=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3; S=double(S);
(3) 复化Cotes公式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); end
for 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; end
C=C+7*(limit(f,a)+limit(f,b)); C=C*h/90; C=double(C);
再编写主程序调用这三个函数,主程序名为ex8_3.m,源程序如下:
%ex8_3.m clc; 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); end
%复化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('三种复化公式积分效果对比图');
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.946083831311699 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