复化求积公式求数值积分.pdf

《复化求积公式求数值积分.pdf》由会员分享，可在线阅读，更多相关《复化求积公式求数值积分.pdf（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、应用软件课程设计 复化求积公式求数值积分 数学 0801 班*、实验目的 程序 1:利用复化梯形公式以及复化辛普森公式求解定积分的数值解。程序 2:分析剖分区间个数对复化梯形公式精度的影响。程序 3:比较 MATLAB 勺 quad、quadl 命令与上述两种方法的精度；比较在相同剖分区间下两种 求数值积分方法的精度。分析与探讨两种方法精度不同的原因。说明：原题目给出的积分精确值 1=4.006994 稍过粗糙，所以通过计算机求解得到更为精确的解 析解。详见测试结果。、算法说明 自定义函数有：2 J-积分函数(在本题中为 I exdx):hanshu.m(附录 1)复化梯形公式求定积分数值解

2、tixing.m(附录 2)复化辛普森公式求定积分数值解 xps.m(附录 3)积分原函数解析解 jxj.m(附录 4)程序 1:p1.m(附录 5)1、n=2000;%确定积分区间分割份数 2、tx=tixing(0,2,n);%用复化梯形公式求解 3、xps=xps(0,2,n);%用复化辛普森公式求解 4、显示结果 disp(积分区间分割分数为：,nu m2str(n)disp(复化梯形公式的求解结果：),disp(tx)disp(复化辛普森公式的求解结果：),disp(xps)程序 2-1:p21.m(附录 6)1、=jxj(2)-jxj(0);%求出解析解。2、tx=zeros(5,

3、1);%给数值解向量赋值。d=zeros(5,1);%给误差向量赋值。3、i=1,2,3,4,5 4、n=109;%定义剖分份数。5、tx(i)=tixing(0,2,n);%将剖分份数 n 代入，求出该 n 下的数值解。d(i)=tx(i)-jxj;%求出误差。6、结束循环 7、输出结果 程序 2-2:p22.m(附录 7)1、=jxj(2)-jxj(0);%求出解析解。2、n=0;l=2000;%迭代次数。dd=5;%迭代步长。d=zeros(1,l);3、i=1,2,1;4、n=n+dd;%积分区间分割份数。d(i)=tix in g(0,2,n)-jxj;%求出误差。5、结束循环 6、

4、figure,x=dd:dd:l*dd;%构造横坐标向量。7、plot(x,d,-m,linewidth,3),xlabel(积分区间等分份数),ylabel(误差大小),title(误差曲线),grid on&%axis(200 10000 0 5*10A-7)%用以控制横纵坐标范围，便于考察不同范围内的误差曲线 情况，可以置于命令窗口中运行。程序 3：p3.m(附录 8)1、a=0;%给定定积分下界。b=2;%给定定积分上界。1=5;%给定数值积分个数。c=1:l;%给定 10 的次方数。n=10.Ac;%给定积分区间的等分份数。jxj=jxj(b)-jxj(a);%求该积分的解析解。tx

5、=zeros(l,1);xp=zeros(l,1);wucha=zeros(l,2);%分别定义两种求数值解的方法不同迭代 次数下的解矩阵以及误差矩阵。2、qwc=abs(quad(hanshu,0,2)-jxj);%求 quad 方法的误差。qlwc=abs(quadl(hanshu,0,2)-jxj);%求 quadl 方法的误差。3、i=c%进行循环 4、tx(i)=tixing(a,b,n(i);%以复化梯形公式求数值解。xp(i)=xps(a,b,n(i);%以复化辛普森公式求数值解。wucha(i,1)=abs(tx(i)-jxj);%求复化梯形公式的误差绝对值。wucha(i,2

6、)=abs(xp(i)-jxj);%求复化辛普森公式的误差绝对值。5、结束循环 6、wc=num2str(c wucha);7、显示运行结果 disp(wc)disp(说明：误差矩阵第一列代表迭代次数 10An。)disp(第二列代表复化梯形公式求得数值解的误差绝对值。)disp(第三列代表复化辛普森公式求得数值解的误差绝对值。)disp(quad 方法求得的解析解误差绝对值为：,n um2str(qwc)disp(quadl 方法求得的解析解误差绝对值为：,n um2str(qlwc)三、测试结果 2 -计算I 0、.一1 exdx的解析解。利用 MATLAB 软件对原函数求不定积分，得到原

7、函数符号解 2 1 ex ln(.ex 1 1)ln(1 1 ex)。在 MATLAB 中 format long 的格式下，显示 I 的值为:4.006994223254704(可在 MATLAB直接通过命令“jxj(2-jxj(0)”驱动 jxj.m 显示该值。)，精 度高于题目所给。为保证计算中保持该数值的最大精度，在求解全过程中，凡是需要此解析解进行 运算的程序，均通过直接调用函数 jxj.m 并将计算结果赋值给变量。运行程序 1,首先将积分区间分割份数 n 取 2000,得到如下结果：积分区间分割分数为：2000 复化梯形公式的求解结果：4.006994300088967 复化辛普森

8、公式的求解结果:4.006994223254719 修改 n 的值，使之等于 10000,得到如下结果:积分区间分割分数为：10000 复化梯形公式的求解结果：4.006994226328080 复化辛普森公式的求解结果：4.006994223254698 运行程序 2-1，得到如下结果:解析解为:4.006994223254704 不同积分区间分割份数下解的值及精度 其中，矩阵第一列分别代表着积分区间分割份数，分割份数分别为：10、102、103、104、105。第二列代表在对应分割份数下使用复化梯形方法得到数值解的值，与该程序输出的解 析解做比较，误差的绝对值显示在矩阵的第三列。由程序 2

9、 运行结果可以得出结论：对于该定积分运算，在复化梯形方法下求解数值解，对积分区间分割份数越多，得到的解越精确。利用 axis 命令可以考察自变量大于 2000 以后的误差随积分区间等分份数之间的关系 运行axis(2000 10000 0 5*10A-8)，得到的函数图像如图 2：1.000000000000000 2.000000000000000 3.000000000000000 4.000000000000000 5.000000000000000 4.010067198653674 4.007024956918338 4.006994530591733 4.0069942263280

10、80 4.006994223285458 0.003072975398970 0.000030733663634 0.000000307337029 0.000000003073376 0.000000000030754 运行程序 2-2，以积分区间等分份数为横轴、数值积分误差为纵轴得到曲线如图 1图1 综合以上结果可以看出，对于该定积分运算，在复化梯形方法下求解数值解，随着积 分区间分割份数的增多，该数值解与解析解的误差呈递减的速率下降，收敛于 0。运行程序 3,得到结果如下：1 0.003073 9.8716e-008 2 3.0734e-005 9.8792e-012 3 3.0734e

11、-007 5.3291e-015 4 3.0734e-009 6.2172e-015 5 3.0754e-011 1.8652e-014 说明：误差矩阵第一列代表迭代次数 10An。第二列代表复化梯形公式求得数值解的误差绝对值。第三列代表复化辛普森公式求得数值解的误差绝对值。quad 方法求得的解析解误差绝对值为：8.9041e-010 quadl 方法求得的解析解误差绝对值为：2.7701e-011 矩阵的第一列意义同程序 2-1，第一列的数字代表 10 的次方数。从该运行结果可以 得出结论：对于该定积分运算，在积分区间等分份数相同的情况下，用复化辛普森公式求 得这个积分的数值解误差绝对值明

12、显小于用复化梯形公式所求的数值解。可以认为，复化 辛普森公式方法的精度高于复化梯形公式法。从输出结果来看，MATLAB的 quadl 命令的计算误差要小于 quad。在积分区间等分 为100 份时，使用复化辛普森公式求得的积分数值解精度已经远高于 quad 和 quadl 方法;在积分区间等分为 100000 份时，使用复化梯形方法下进行运算时其精度高于 quad 方法，虽接近但并未高于 quadl 算法的精度。图2 四、分析与探讨 几何上，图 3 中辛普森方法下多面体面积（黑框内部分）相比梯形方法下更加接近积分面积真实 值，因此相同积分区间相同分割份数下，前者的计算结果要更加准确。当积分区间

13、等分份数超过某个临界值时，两种算法精度可以分别优于 quad、quadl 算法。根 据以上分析，复化辛普森方法的此临界值应当小于复化梯形方法的。为了得到该临界值并验证上 述结论的正确性，现编写了程序 analyse.m（参见附录 9,因原理相同，四个程序的源代码仅列举 一个的）。运行 analyse，通过不断修改初值、迭代，最终得到结果如下：使用 quad 数值积分命令的误差为：8.9041e-010 等分份数：18578 误差：8.9045e-010 等分份数：18579 误差：8.9038e-010 为使得复化梯形方法精度高于 quad，区间均分份数至少为：18579 使用 quadl 数

14、值积分命令的误差为：2.7701e-011 等分份数：105269 误差：2.7736e-011 由题目知，复化梯形方法公式:b f(x)dx a 1 丄(b a)f(a)f(b)，2 复化辛普森公式:b b a a b af(x)dx f(a)4f(亍)f(b)1 a b-f(a)4f(-)f(b)比复化梯形公式在改点取值 6 2 由于复化辛普森公式在 咛更加接近 x-点的取值 2 f(-b)(如图 3)。在 2 等分份数：105270 误差：2.7689e-011 为使得复化梯形方法精度高于 quadl，区间均分份数至少为：105270 使用 quad 数值积分命令的误差为：8.9041e

15、-010 等分份数：32 误差：9.4199e-010 等分份数：33 误差：8.329e-010 为使得复化辛普森方法精度高于 quad，区间均分份数至少为：33 使用 quadl 数值积分命令的误差为：2.7701e-011 等分份数：77 误差：2.8103e-011 等分份数：78 误差：2.6687e-011 为使得复化辛普森方法精度高于 quadl，区间均分份数至少为：78 从上述程序运行结果很容易验证上述结论的正确性。相比复化梯形公式，使用复化辛普森公式 求解定积分具有无可比拟的优越性。不可否认，相同分割区间份数下复化辛普森方法要比复化梯形方法复杂。然而，如果放开分割 分数相同这

16、一条件，在精度要求相同时前者运行时间应当小于后者。为了验证这一结论，现编写了 程序text.m（参见附录 10），考察两种复化方法在精度分别高于 quad 方法和 quadl 方法时程序的 运行速度。程序运行结果如下：quad 方法 Elapsed time is 0.006125 sec on ds.quadl 方法 Elapsed time is 0.048047 seco nds.复化梯形方法迭代 18590 步（与 quad 比较）Elapsed time is 0.843971 seco nds.复化梯形方法迭代 105370 步（与 quadl 比较）Elapsed time is

17、 3.515906 sec on ds.复化辛普森方法迭代 33 步：（与 quad 比较）Elapsed time is 0.001900 seco nds.复化辛普森方法迭代 78 步：（与 quadl 比较）Elapsed time is 0.003362 sec on ds.从表 1 中不难看出，在复化梯形方法、复化辛普森方法分别与qual、quadl 两种方法精度相同 时，复化辛普森方法的运行时间都远小于复化梯形方法的。于是前一结论得到验证。有关误差的解释：有同学提出，原积分是符号解，在 MATLAB 总是存在误差的，因此对建立 在该误差上的 quad、quadl 精度讨论无意义。但

18、本文给出的积分精确值保留了 15 位小数（误差水 方法 quad quadl 复化梯形方法 复化辛普森方法 比较项目 误差 运行时间 误差 运行时间 误差 运行时间 误差 运行时间 值 89.041 0.006125 89.038 0.843971 83.29 0.001900 值 2.7701e 0.048047 2.7689 3.515906 2.6687 0.003362 将以上 analyse.m、text.m 程序运行结果汇总到如下表 1 中:注：误差的单位为 10 11 平 1 10 15），而对 quad、quadl 的讨论误差在1 10 10 1 10 11水平上，因此认为 M

19、ATLAB!身因 符号解转化为数值产生的误差对该讨论不产生影响，因此忽略不计。附录：1、hanshu.m function y=hanshu(x)y=sqrt(1+exp(x);2、tixing.m function s=tixing(a,b,n)%a:下界，b:上界，n：区间分割份数。ss=0;d=(b-a)/n;x=a:d:b;for i=1:n-1;ss=ss+hanshu(a+i*d);%end fa=hanshu(a);fb=hanshu(b);s=d*(1/2)*(fa+2*ss+fb);%加和求复化梯形面积。End 3、xps.m function s=xps(a,b,n)%a:

20、下界，b:上界，n：迭代次数。s=hanshu(a)+hanshu(b);h=(b-a)/n;for i=1:n-1;s=s+2*hanshu(a+i*h);end for k=1:n;s=s+4*hanshu(a+h*(k-0.5);end s=s*h/6;4、jxj.m function jxj=jxj(x)%输入自变量 x 的值，输出该问题原函数的解析解。u=sqrt(1+exp(x);jxj=2.*u+log(u-1)-log(1+u);5、p1.m clc,clear,format longn=2000;%积分区间分割份数 tx=tixing(0,2,n);%用复化梯形公式求解 xp

21、s=xps(0,2,n);%用复化辛普森公式求解 disp(积 分 区 间 分 割 分 数 为:,num2str(n),disp(复 化 梯 形 公 式 的 求 解 结 果:),disp(tx),disp(复化辛普森公式的求解结果:),disp(xps)6、p21.m clc,clear,format long jxj=jxj(2)-jxj(0);%求出解析解。tx=zeros(5,1);%给数值解向量赋值。d=zeros(5,1);%给误差向量赋值。for i=1:5;n=10Ai;%定义剖分份数。tx(i)=tixing(0,2,n);%将剖分份数 n 代入，求出该 n 下的数值解。d(i

22、)=tx(i)-jxj;%求出误差。end j=1:5;disp(解析解为：)disp(jxj),disp(不同积分区间分割份数下解的值及精度),disp(j tx d);7、p22.m clc,clear,format long jxj=jxj(2)-jxj(0);%求出解析解。n=0;d=zeros(1,l);1=2000;%迭代次数。dd=5;%迭代步长。for i=1:l;n=n+dd;%积分区间分割份数。d(i)=tixing(0,2,n)-jxj;%求出误差。end figure x=dd:dd:l*dd;%构造横坐标向量。plot(x,d,-m,linewidth,3),xlab

23、el(积分区间等分份数),ylabel(误差大小),title(误差 曲线),grid on%axis(200 10000 0 5*10A-7)%用以控制横纵坐标范围，便于考察不同范围内的误差曲线情况，可 以置于命令窗口中运行。8、p3.m clc,clear format long%将输出显示格式调整到小数点后 14 位，便于对比微小误差。a=0;%给定定积分下界。b=2;%给定定积分上界。1=5;%给定数值积分个数。c=1:l;n=10.Ac;jxj=jxj(b)-jxj(a);%求该积分的解析解。tx=zeros(l,1);xp=zeros(l,1);wucha=zeros(l,2);%

24、分别定义两种求数值解的方法不同迭代次数下 的解矩阵以及误差矩阵。qwc=abs(quad(hanshu,0,2)-jxj);%求 quad 方法的误差。qlwc=abs(quadl(hanshu,0,2)-jxj);%求 quadl 方法的误差。for i=c tx(i)=tixing(a,b,n(i);%以复化梯形公式求数值解。xp(i)=xps(a,b,n(i);%以复化辛普森公式求数值解。wucha(i,1)=abs(tx(i)-jxj);%构造误差矩阵中复化梯形公式的误差绝对值。wucha(i,2)=abs(xp(i)-jxj);%构造误差矩阵中复化辛普森公式的误差绝对值。end wc

25、=num2str(c wucha);disp(wc),disp(说明：误差矩阵第一列代表迭代次数 10An。)disp(第二列代表复化梯形公式求得数值解的误差绝对值。)disp(第三列代表复化辛普森公式求得数值解的误差绝对值。)disp(quad 方法求得的解析解误差绝对值为：,num2str(qwc)disp(quadl 方法求得的解析解误差绝对值为：,num2str(qlwc)9、analyse.m(因原理相同，四个程序的源代码仅列举第一个)clc,clear;format long%将输出显示格式调整到小数点后 14 位，便于对比微小误差。a=0;%给定定积分下界。b=2;%给定定积分上

26、界。dd=1;%迭代步长。ll=18570;%给定数值积分区间等分份数下界。lu=18590;%给定数值积分区间等分份数上界。n=ll:dd:lu;jxj=jxj(b)-jxj(a);%求该积分的解析解。qwc=abs(quad(hanshu,0,2)-jxj);disp(使用 quad 数值积分命令的误差为：,num2str(qwc);for i=1:length(n);tx=tixing(a,b,n(i);%以复化梯形公式求数值解。aa=abs(tx-jxj);%求复化梯形公式的误差。disp(等分份数：,num2str(n(i),误差：,num2str(aa)if aa=qwc if i

27、=1;disp(等分份数已足够大，在下限处复化梯形方法精度就已高于 quad！)return elseif ilength(n);disp(为使得复化梯形方法精度高于 quad,区间均分份数至少为：,num2str(n(i)return elseif i=length(n);disp(在该等分份数下,复化梯形方法精度始终低于 quad！)end end end 10、text.m xps(0,2,78);toc clc,clear disp(quad 方法)tic quad(hanshu,0,2);toc disp(),tic quadl(hanshu,0,2);toc disp(),tic tixing(0,2,18590);toc disp(),tic tixing(0,2,105370);toc disp(),tic xps(0,2,33);toc disp(),tic disp(quadl 方法)disp(复化梯形方法迭代 18590 步(与 quad 比较)disp(复化梯形方法迭代 105370 步(与 quadl 比较)disp(复化辛普森方法迭代 33 步：(与 quad 比较)disp(复化辛普森方法迭代 78 步：(与 quadl 比较)