复化梯形复化辛普森龙贝格自适应辛普森课程设计论文

湖南农业大学综合设计报告综合设计五多方法求解数值积分学生姓名：学号：年级专业：指导老师：学院：评阅成绩：评阅意见：成绩评定教师签名：时间：湖南·长沙提交日期：2014年6月多方法求解数值积分具体题目要求：用不同数值方法计算积分104ln9xxdx(1)取不同的步长h，分别用复合梯形及复合辛普森公式计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？(2)用龙贝格求积计算完成问题(1);(3)用自适应辛普森积分，使其精度达到410。1设计目的、要求由积分学基本理论，定积分可由NewtonLeibniz公式计算，但是对于一些无法找到原函数的函数（如2xe等）不能通过牛顿—莱布尼兹公式计算，就必须得另寻它法。因此需要我们能够熟练地应用常用的数值积分计算方法（如机械求积、NewtonCotes公式等）并掌握结合数值计算软件（Matlab、Lingo等）及计算机高级语言()cjava、进行对应算法实现的技能。熟练数学软件求解数学问题，掌握各种数学问题的求解方法。本设计主要是通过多种复合求积公式求解积分，主要包括复化梯度法、复化辛普森法、龙贝格以及自适应辛普森法等求解方法，利用Matlab软件编写相对应的算法进行求解，大大地提高了解题的速度。2设计原理由积分中值定理我们可以知道在积分区间,ab内存在一点，使得式子()()()bafxdxbaf成立。这个式子在于对于点的具体位置一般是不知道的，因此难以准确算出()f的值。也就是不同算法求得平均高度()f，对应的就是一种不同的数值求积方法。更一般地，我们可以在区间,ab上适当选取某些节点kx，然后用的加权平均得到平均高度()f的近似值，这样构造出的求积公式具有下列形式：0()()nbkkakfxdxAfx称为机械求积公式。复合梯形公式、复合辛普森公式、龙贝格求积公式以及自适应辛普森公式都以此公式的基础，对积分区间进行变步长的划分求得近似的平均高度值，得到积分函数的近似值。也由于牛顿—柯特斯公式在8n时不具有稳定性，所以不可能通过提高阶的方法来提高求积精度。为了我提高精度通常可把积分区间等分成为若干个子区间，再在每个子区间上用低价求积公式，这就是复合求积方法。但是这样的积分求解方法也是存在不容忽视的误差。因此需要在设计算法时考虑到算法存在的误差（舍入误差、截断误差等），并对误差作出分析。3采用软件、设备Matlab软件4设计内容第一步：复合梯形公式、复合辛普森公式算法（一）、复合梯形公式计算积分复化梯形公式的主要思想是利用若干小梯形的面积代替原方程的积分，利用微元法，可以求出坐标面上由函数与坐标轴围城的图像的面积的近似值，符合了计算机计算存储的思想。下面，我们在探讨复化梯形公式的计算规律:设将求积区间ba,分成n等份，则一共有1n个分点，按梯形公式bfxfafhxfxfhTnkkkknkn11110222计算积分值nT，需要提供1n个函数值。这里nabh代表步长，分点为khaxk，其中.,1,0nk（二）、复合辛普森公式计算积分算法的基本思想是：把积分区间等分成若干个子区间，而在每一个子区间上用辛普森求积公式：[()4()()]62baabSfaffb得到复合辛普森求积公式：111012[()4()2()()]6nnnkkkkhsfafXfXfb并且用Matlab软件来求解。第二步：龙贝格算法考虑积分()()baIffxdx，欲求其近似值，通常有复化的梯形公式、Simpsion公式和Cotes公式。但是给定一个精度，这些公式达到要求的速度很缓慢。如何提高收敛速度，自然是人们极为关心的课题。为此，记1T，k为将区间,ab进行2k等分的复化的梯形公式计算结果，记2T，k为将区间,ab进行2k等分的复化的Simpsion公式计算结果，记3T,k为将区间,ab进行2k等分的复化的Cotes公式计算结果。根据Richardson外推加速方法，可以得到收敛速度较快的Romberg积分法。其具体的计算公式为：2,mk，1、准备初值，计算)]()([21,1bfafbaT2、按梯形公式的递推关系，计算1201,11,11))5.0(2(221kikkkkiabafabTT3、按Romberg积分公式计算加速值11,11,,1441mmkmmkmmkmmTTT，2,mk，4、精度控制。对给定的精度R，若RTTmm1,11,则终止计算，并取,1mT为所求结果；否则返回2重复计算，直至满足要求的精度为止。第三步：自适应辛普森算法复合求积方法通常适用于被积函数变化不太大的积分，如果在积分区间被积函数变化很大有的部分函数值变化剧烈而有的部分则是变化平缓，如果此时还是将积分区间等分后用复合求积公式的话计算量很大。而采用针对被积函数在区间上的不同情形而采用不同的步长，使得在满足精度的前提下积分计算量减少，这就是自适应积分方法，能自动地在被积函数变化剧烈的区域增多节点，而在被积函数变化平缓的地方减少节点。因此它是一种不均匀区间的积分方法。题目要求使相邻两个区间的误差达到一定的要求410，即用自适应辛普森公式来求积分，先算出积分区间的左右端点函数值，求出区间中点函数值与左右端点的函数差值，再与所要求的精度比较，不满足的对所在区间二等分，接着算出每个子区间端点的函数值判断时否符合精度要求，直到积分每个子区间内都满足精度要求，最后所得各个区间端点的函数值之和即为积分的近似值。第四步：误差余项以及精度分析：由插值型的求积公式我们得到求积公式误差余项的表达式：10nbmkkakRffxdxAfxKf其中m表示求积公式的代数精度，为不依赖于fx的待定系数，(,)ab，这个结果表明当()fx是次数小于等于m的多项式时由于(1)()0mfx，故此时[]0Rf，即前面的求积公式精确成立。而当1()mfxx时，(1)()(1)!mfxm故可求得：1122100111112nnbmmmmmkkkkakkKxdxAxbaAxmmm因为梯形公式的代数精度为1，可以得到K的值为33322111123226baKbaabba于是得到梯形公式的余项为：3......,12baRffab又因为复合梯形公式要满足101nkkffn综上所述，就得到了复合梯形公式的余项表达式：2......,12nbaRfhfab同理可得复合辛普森公式的余项表达式：41410......,1802nnnkkkkkbahRfISfxx结果分析：从以上余项的表达式可以看出复合辛普森公式的代数精度为3，而复合梯形公式的代数精度为1，所以复合辛普森比复合梯形精确度更高。对于算法的精度，是通过对设计所得值与准确值之间的误差值来评判。将变步长的复合求积方法每次求得计算结果与准确值进行比较求出误差值，通过画出误差值的变化趋势图比较复合梯形公式与复合辛普森公式这两种算法的精度。经过实验中验证，也表明自己的初步推理是正确的，无论是复合梯形公式还是复合辛普森公式它们最终结果都会随着步长值的减小而更加精确。复合梯形公式和复合辛普森公式计算出的结果进行比较，发现复合辛普森公式计算出的结果更加的精确。5原始程序、数据文件：f.mfunctiony=f(x);y=sqrt(x)*log(x);1、复合梯形公式求解算法：文件：trapezoid.mclca=0;%积分下限b=1;%积分上限T=[];%用来装不同n值所计算出的结果R=[];G=[];m=120;%等分数true=-(4/9);forn=2:m;h=(b-a)/n;%步长x=zeros(1,n+1);%给节点定初值y=zeros(1,n+1);fori=1:n+1x(i)=a+(i-1)*h;endx(1,1)=0.000000001;fori=1:n+1y(i)=x(i).^(1/2)*log(x(i))%g=(-(b-a)/12*h^2)*(-log(x(i))/(4*x(i)*x(i)^(1/2)))%准确的积分余项（计算误差）end%G=[G,g];t=0;r=0;fori=1:nformatlongt=t+h/2*(y(i)+y(i+1));%利用复化梯形公式求值err=t-floor(t);digits(7);%此处为需要的小数位+1t=floor(t)+vpa(err,6);%此处控制显示的小数点位数,更改显示的小数位数r=t-true;%计算的值与真实值之差（实际误差）endT=[T,t];%把不同n值所计算出的结果装入T中R=[R,r];endx=linspace(0,1,m-1);plot(x,R,'*')%将计算误差与实际误差用图像画出来2、复合辛普森积分求解算法：simpon.mclccleara=0;%积分下限b=1;%积分上限T=[];%用来装不同n值所计算出的结果R=[];true=-(4/9);m=20;%等分数forn=2:mh=(b-a)/(2*n);%步长x=zeros(1,2*n+1);%给节点定初值y=zeros(1,2*n+1);fori=1:2*n+1x(i)=a+(i-1)*h;%给节点赋值endx(1,1)=0.000000001;fori=1:2*n+1y(i)=x(i).^(1/2)*log(x(i));%给相应节点处的函数值赋值endt=0;r=0;fori=1:nformatlongt=t+h/3*(y(2*i-1)+4*y(2*i)+y(2*i+1));%利用复化simpson公式求值err=t-floor(t)digits(7);%此处为需要的小数位+1t=floor(t)+vpa(err,6);r=t-true;endT=[T,t]%把不同n值所计算出的结果装入T中R=[T,r];end%R=(-(b-a)/180*((b-a)/2).^4*24)%积分余项（计算误差）%true=quad(@fx1,0,1);%积分的真实值x=linspace(0,1,2*m-1);plot(x,R,'*')3、龙贝格算法rebeg.m%%龙贝格clearclca=0;b=1;%确定积分上下限eps=10^(-4);err=1;k=1;a=0.0000001;T(1,1)=(b-a)/2*(f(a)+f(b));while(erreps)h=(b-a)/2^(k-1);S=0;forx=a:h:b-hS=S+f(x+h/2);endT(k+1,1)=1/2*T(k,1)+h/2*S;k=k+1;fori=2:kT(k,i)=(4^(k-1)*T(k,i-1)+T(k-1,i-1))/(4^(k-1)-1);enderr=abs(T(i,i)+4/9);endfprintf('龙贝格求积算法积分值为%.10f\n',T(k,k));disp(T)T龙贝格求积算法积分值为-0.44438207534、自适应辛普森算法：%自适应辛普森算法Self_Adaptive_integral.mfunctions=Self_Adaptive_integral(a,b,tol)k=0;w=0;x=a;y=b;t=0;h=(b-a)/2;s=0;i=0;to=abs(simpson_integral(x,y,2)-simpson_integral(x,y,1));whileto=toli=i+1;whileto=tolt=x;ifk==0x=t;y=t+h;to=(abs