0.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.946083071968057 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.946020325355025 0.946083070751532 0.946083070367182
三种复化公式积分效果对比图0.955准确值曲线复化梯形公式复化Simpson公式复化Cotes公式 0.950.9450.940.9350.930.9250.92 02468101214161820 由以上结果可看到复化梯形公式有一个上升接近准确值过程,而复化Simpson公式积分结果和复化Cotes公式积分的结果基本上和准确值的曲线重叠在一块,可见它们的精度是相当高的.
例 8.4 用Romberg积分法计算I??解：
编写Romberg积分法的函数M文件,源程序如下(romberg.m):
function [I,T]=romberg(f,a,b,n,Eps) %Romberg积分计算 %f为积分函数 %[a,b]为积分区间
sinxdx,精度??10?6. 0x1%n+1是T数表的列数目 %Eps为迭代精度
%返回值中I为积分结果,T是积分表
if nargin<5 Eps=1E-6; end m=1; h=(b-a); err=1; j=0;
T=zeros(4,4);
T(1,1)=h*(limit(f,a)+limit(f,b))/2; while ((err>Eps) & (j for p=1:m x0=a+h*(2*p-1); s=s+limit(f,x0); end T(j+1,1)=T(j,1)/2+h*s; m=2*m; for k=1:j T(j+1,k+1)=T(j+1,k)+(T(j+1,k)-T(j,k))/(4^k-1); end err=abs(T(j,j)-T(j+1,k+1)); end I=T(j+1,j+1); if nargout==1 T=[]; end 将上述源程序另存为romberg.m后,进入计算: >>syms x;%创建符号变量 >>f=sym('sin(x)/x') %符号函数 f = sin(x)/x >> [I,T]=romberg(f,0,1,3,1E-6) %积分计算I = 0.9461 T = 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为Romberg积分表,由输出结果可知计算结果为I?sinx?0xdx?0.9461. 1 例 8.5 利用Romberg积分法求I??1401?x2dx. 解: 直接利用例8.4的函数即可. >>syms x; >>f=sym('4/(1+x^2)') f = 4/(1+x^2) >>[I,T]=romberg(f,0,1,5,1E-6) I = 3.1416 T = 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 故结果为I??1401?x2dx?3.1416. 例 8.6 对积分 ?1?1f(x)dx?A0f(x0)?A1f(x1) 构造其Gauss型求积公式. 解: 取f(x)?1,x,x2,x3,得到方程组: ??A0?A1?2??A0x0?A1x1?0??A0x20?A1x21?2 ?3??A330x0?A1x1?0 这是一个抽象方程组,可以利用Matlab的符号法来解之,该函数名为solve(): %直接输入解之: >>x=solve('A0+A1=2','A0*x0+A1*x1=0','A0*x0^2+A1*x1^2=2/3','A0*x0^3+A1*x1^3=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 1 ans = 1 1 ans = 1/3*3^(1/2) -1/3*3^(1/2) ans = -1/3*3^(1/2) 1/3*3^(1/2) 因此得到两组解为: (A0,A1,x0,x1)?(1,1,33,?) 3333,) 33或: (A0,A1,x0,x1)?(1,1,?求积公式为 ? 1?1f(x)dx?f(33)?f(?) 33例 8.13 用中心差商公式计算f(x)?解： 应用中心差商公式可得到： x在x?2处的一阶导数. 2?f'(2)?编制如下程序计算： hh?2?22 h %ex8_13.m clc; syms x;%符号变量 f=sym('sqrt(x)')%符号函数 x0=2;%求值点 h=2;%第一次初始步长 N=10;%进行10次迭代计算 dF=diff(f,x,1)%求导 dFx=subs(dF,x0)%精确值 df=zeros(N,6); for i=1:N df(i,1)=h; temp=(subs(f,x,x0+h/2)-subs(f,x,x0-h/2))/h;%按中心差商公式计算 df(i,2)=double(temp); h=h/10;%步长缩小10倍 end temp=dFx*ones(N,1)-df(:,2); df(:,3)=temp; h=1;%第二次取初始步长为1 for i=1:N df(i,4)=h; temp=(subs(f,x,x0+h/2)-subs(f,x,x0-h/2))/h; df(i,5)=double(temp); h=h/10; end temp=dFx*ones(N,1)-df(:,5); df(:,6)=temp; disp(' 步长h 近似值 误差 步长h 近似值 误差 '); disp(df); 运行结果如下： dFx = 0.353553390593274 步长h 近似值 误差 步长h 近似值 误差 Columns 1 through 4 2.000000000000000 0.366025403784439 -0.012472013191165 1.000000000000000 0.200000000000000 0.353663997049609 -0.000110606456335 0.100000000000000 0.020000000000000 0.353554495459696 -0.000001104866422 0.010000000000000 0.002000000000000 0.353553401641782 -0.000000011048508 0.001000000000000 0.000200000000000 0.353553390703976 -0.000000000110702 0.000100000000000