8数值积分与数值微分

8数值积分与数值微分8.1例题解答例8.1给定积分10.5xdx,分别用梯形公式、Simpson公式、Cote公式作近似计算.解:先输入主要初始参数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=[73212327];fori=0:4tc=tc+C0(i+1)*feval(f,a+i*(b-a)/4);endI3=(b-a)/90*tcI3=0.430964070495876%准确值I=int(char(f),a,b)vpa(I)I=-1/6*2^(1/2)+2/3ans=0.43096440627115082519971854596505例8.2对积分10sinxIdxx,为使其精度达到410.若用复化梯形公式,应将[0,1]多少等分?若用复化Simpson公式,应将[0,1]多少等分?解:直接按余项计算即可.复化梯形公式的余项为:2''()(),(,)12nbaEfhfab复化Simpson公式余项为:4(4)4()(),(,),()[,]1802nbahEffabfxCab对于sin()xfxx,在课本中我们已证得以下不等式成立:()1()1kfxk直接利用上述不等式关系解答本题.先输入误差精度:Eps=1E-4Eps=1.000000000000000e-004(1)复化梯形公式h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1)))%先求出步长h1=0.060000000000000N1=ceil(1/h1)%向上取整,得到等分区间数N1=17故可将区间17等分即可达到所要求的精度.(2)复化Simpson公式h2=power(Eps/abs(-(1-0)/180*1/(1+4)),1/4)%先求出步长h2=0.547722557505166N2=ceil(1/h2)%向上取整,得到等分区间数N2=2故可将区间2等分即可达到所要求的精度.扩展:1)Matlab中复化梯形公式命令为I=trapz(x,y),复化Simpson公式命令为quad().2)Matlab中有四个取整函数,分别为ceil(),floor(),fix(),round(),分别表示向正无穷大方向取整、向负无穷大方向取整、向靠近零方向舍入和四舍五入.例8.3对积分10sinxIdxx,利用变步长方法求其近似值,使其精度达到610.解：利用变步长法前先建立三种变步长复化积分公式的函数.注意在Matlab中直接用sin(0)/0得不到1，sin(0)0matlabNAN，因此解此题时我们改用求极限的方法得到函数值,此函数名为limit().先建立三种复化公式的函数文件,它们分别为复化梯形公式trap.m、复化Simpson公式为simpson.m、Cotes公式为cotes.m,三个函数的源程序如下:(1)复化梯形公式trap.mfunctionT=trap(f,a,b,n)%trap.m%复化梯形公式求积分值%f为积分函数%[a,b]为积分区间%n是等分区间份数h=(b-a)/n;%步长T=0;fork=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公式simpson.m:functionS=simpson(f,a,b,n)%simpson.m%Simpson公式求积分值%f为积分函数%[a,b]为积分区间%n是等分区间份数h=(b-a)/(2*n);%步长s1=0;s2=0;fork=1:nx0=a+h*(2*k-1);s1=s1+limit(f,x0);endfork=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公式cotes.m:functionC=cote(f,a,b,n)%cote.m%复化cotes公式求积分值%f为积分函数%[a,b]为积分区间%n是等分区间份数h=(b-a)/n;%步长C=0;fori=1:(n-1)x0=a+i*h;C=C+14*limit(f,x0);endfork=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;symsx;f=sym('sin(x)/x');a=0;b=1;%积分上下限n=20;%作1,2,3,…,20次区间等分%复化梯形公式T=zeros(n,1);fori=1:nT(i)=trap(f,a,b,i);end%复化Simpson公式;S=zeros(n,1);fori=1:nS(i)=simpson(f,a,b,i);end%复化Cotes公式C=zeros(n,1);fori=1:nC(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,'-');holdon;plot(x,T','r--','LineWidth',2);plot(x,S','m.-','LineWidth',1);plot(x,C','c:','LineWidth',1.5);gridon;title('三种复化公式积分效果对比图');legend('准确值曲线','复化梯形公式','复化Simpson公式','复化Cotes公式');holdoff;disp('复化梯形公式复化Simpson公式复化Cotes公式');disp([TSC]);在Matlab命令窗口中输入ex8_3,得到如下积分结果:复化梯形公式复化Simpson公式复化Cotes公式0.9207354924039480.9461458822735870.9460830040636740.9397932848061770.9460869339517940.9460830693509170.9432914291323370.9460838313116990.9460830702782780.9445135216653900.9460833108884720.9460830703513790.9450787809534020.9460831688380730.9460830703630430.9453857307668590.9460831178428670.9460830703657970.9455707762562460.9460830959894030.9460830703666330.9456908635827010.9460830853849480.9460830703669360.9457731885497520.9460830797420530.9460830703670610.9458320718669050.9460830765177320.9460830703671180.9458756371191980.9460830745679370.9460830703671470.9459087710738650.9460830733331140.9460830703671610.9459345564884750.9460830725204760.9460830703671700.9459550160561140.9460830719680570.9460830703671740.9459715215529850.9460830715819640.9460830703671770.9459850299343860.9460830713055620.9460830703671790.9459962252423760.9460830711034890.9460830703671800.9460056069439780.9460830709529980.9460830703671810.9460135466239660.9460830708390650.9460830703671820.9460203253550250.9460830707515320.946083070367182024681012141618200.920.9250.930.9350.940.9450.950.955三种复化公式积分效果对比图准确值曲线复化梯形公式复化Simpson公式复化Cotes公式由以上结果可看到复化梯形公式有一个上升接近准确值过程,而复化Simpson公式积分结果和复化Cotes公式积分的结果基本上和准确值的曲线重叠在一块,可见它们的精度是相当高的.例8.4用Romberg积分法计算10sinxIdxx,精度610.解：编写Romberg积分法的函数M文件,源程序如下(romberg.m):function[I,T]=romberg(f,a,b,n,Eps)%Romberg积分计算%f为积分函数%[a,b]为积分区间%n+1是T数表的列数目%Eps为迭代精度%返回值中I为积分结果,T是积分表ifnargin5Eps=1E-6;endm=1;h=(b-a);err=1;j=0;T=zeros(4,4);T(1,1)=h*(limit(f,a)+limit(f,b))/2;while((errEps)&(jn))|(j4)j=j+1;h=h/2;s=0;forp=1:mx0=a+h*(2*p-1);s=s+limit(f,x0);endT(j+1,1)=T(j,1)/2+h*s;m=2*m;fork=1:jT(j+1,k+1)=T(j+1,k)+(T(j+1,k)-T(j,k))/(4^k-1);enderr=abs(T(j,j)-T(j+1,k+1));endI=T(j+1,j+1);ifnargout==1T=[];end将上述源程序另存为romberg.m后,进入计算:symsx;%创建符号变量f=sym('sin(x)/x')%符号函数f=sin(x)/x[I,T]=romberg(f,0,1,3,1E-6)%积分计算I=0.9461T=0.920700000.93980.94610000.94450.94610.9461000.94570.94610.94610.946100.94600.94610.94610.94610.9461其中T为Romberg积分表,由输出结果可知计算结果为10sin0.9461xIdxx.例8.5利用Romberg积分法求12041Idxx.解:直接利用例8.4的函数即可.symsx;f=sym('4/(1+x^2)')f=4/(1+x^2)[I,T]=romberg(f,0,1,5,1E-6)I=3.1416T=3.0000000003.10003.133300003.13123.14163.14210003.13903.14163.14163.1416003.14093.14163.14163.14163.141603.14143.1