合肥工业大学计算方法总结计划复化梯形公式实验.docx
- 格式:docx
- 大小:61.61 KB
- 文档页数:9
《计算方法》实验报告(2) 虽然找到了原函数,但因表达式过于复杂而不便计算 (3) f(x)是由测量或计算得到的表格函数由于以上种种困难,有必要研究积分的数值计算问题。
利用插值多项式 )()(x f x P n ≈则积分dx x f ba⎰)(转化为dx x P b a n ⎰)(,显然易算。
dx x P ba n ⎰)(称为插值型求积公式。
最简单的插值型求积公式是梯形公式和Simpson 公式,。
当求积结点提供较多,可以分段使用少结点的梯形公式和Simpson 公式,并称为复化梯形公式、复化Simpson 公式。
如步长未知,可以通过误差限的控制用区间逐次分半的策略自动选取步长的方法称自适应算法。
梯形递推公式给出了区间分半前后的递推关系。
由梯形递推公式求得梯形序列,相邻序列值作线性组合得Simpson 序列, Simpson 序列作线性组合得柯特斯序列, 柯特斯序列作线性组合的龙贝格序列。
若|R 2-R 1|<,则输出R2;否则…依此类推。
如此加工数据的过程叫龙贝格算法,如下图所示:复化梯形公式)]()(2)([21110n n k k n x f x f x f h T ++=∑-=复化Simpson 公式])(4)()(2)([61021110∑∑-=+-=+++=n k k n k n k n x f x f x f x f hS梯形递推公式∑-=++=10212)(22n k k n n x f h T T 加权平均公式:n n n S T T =--1442 n n n C S S =--144222 n nn R C C =--144323 龙贝格算法大大加快了误差收敛的速度,由梯形序列O(h2) 提高到龙贝格序列的O(h8)return R2;}}实验结果:图1六、实验总结1.梯形公式的收敛速度太慢,所以我们才会选择之后的几种公式加快收敛速度。
2.通过实验可以看出,龙贝格算法大大加快了误差收敛的速度,由梯形序列O(h2) 提高到龙贝格序列的O(h8)3.在编程的过程中可以切身体验到龙贝格算法是基于前几种算法的实现才能实现的。
第二章1.1复合梯形求积公式复合梯形求积公式是复合求积法的一种,在本章中,将从其原理、概念等方面对它做一个详细介绍。
在本章的最后,会对复合梯形求积法进行程序设计,使得可以从不同的方面对这种方法有更深的理解。
1.1.1 复合梯形求积公式的理论当积分区间[a ,b]的长度较大,而节点个数1+n 固定时,直接使用Newton-Cotes 公式的余项将会较大。
但是如果增加节点个数,即1+n 增加时,公式的舍入误差又很难得到控制。
为了提高公式的精度,又使算法简单易行,往往使用复化方法。
即将积分区间][b a ,分成若干子区间,然后在每个小区间上使用低阶Newton-Cotes 公,最后将每个小区间上的积分的近似值相加,这就叫做复合求积法。
而复合梯形求积公式就是复合求积法的一种。
1.1.2复合求积公式的原理将区间[]b a ,划分为n 等分,分点,,,1,0,,n k nab h kh a x k =-=+= 在每个子区间[](),1,,1,0,1-=+n k x x k k 上采用梯形公式,则得[])()()(2)()(1111f R x f x f h dx x f dx x f I n k n k k ban k x x k k++===+-=-=∑⎰∑⎰+记()[()]()[()()]∑∑-=+-=++=+=11110222n k b k k n k k n x f x f a f hx f x f h T , (1.1)称为复合梯形公式,其余项可由)().,(),(12][''3b a f a b f R ∈--=ηη得()()()110''3,12+-=∈⎥⎦⎤⎢⎣⎡-=-=∑k k k n k k n n x x f h T I f R ηη由于[],,)(2b a C x f ∈ 且()(),max 1min 1010''''10-≤≤-=-≤≤≤≤∑n k k n k k n k f n f ηη所以()b a ,∈∃η使()()k n k f n f ηη∑-==10''''1于是复合梯形公式的余项为()()η''212f h a b f R n --= 可以看出误差是2h 阶,且由()()η''212f h a b f R n --= 式立即得到,当][b a C x f ,)(2∈时,则(),lim dx x f T ba n n ⎰=∞→即复合梯形公式是收敛的.事实上只要设()[]b a C x f ,∈,则可得收敛性,只要把n T 改写成为()()]∑∑=-=-+⎢⎣⎡-=nk k n k k n x f n a b x f n a b T 11021当∞→n 时,上式右端括号内的两个和式均收敛到积分()dx x f ba ⎰,所以复合梯形公式(1.1)收敛.此外,n T 的求积系数为正,由定理可知复合梯形公式是稳定的。
利用复化梯形公式复化公式计算积分设被积函数为f(x),要在区间[a,b]上进行积分,将区间[a,b]均分为n个小区间,每个小区间的宽度为h=(b-a)/n。
则,复合梯形公式用来计算区间[a,b]上的积分可以表示为:∫[a, b] f(x) dx ≈ h / 2 * [ f(a) + 2f(x1) + 2f(x2) + ... + 2f(xn-1) + f(b) ]其中,x1, x2, ..., xn-1为子区间的起始点,满足 xi = a + i * h,i = 1, 2, ..., n-1复合梯形公式的误差估计为E(h)=-(b-a)*h²/12*f''(ξ),其中ξ∈[a,b],f''(ξ)为f(x)的二阶导数。
下面以一个具体例子来说明如何利用复合梯形公式进行积分计算。
例如,要计算函数f(x)=x²在区间[0,1]上的积分。
首先,选择将区间[0,1]均分为n个小区间,每个小区间的宽度为h=(1-0)/n。
根据复合梯形公式,我们可以得到积分的近似值为:∫[0, 1] x² dx ≈ h / 2 * [ f(0) + 2f(x1) + 2f(x2) + ... + 2f(xn-1) + f(1) ]其中,xi = 0 + i * h,i = 1, 2, ..., n-1这样,我们可以计算出每个子区间的积分值,并将它们加总得到最终的结果。
例如,当n=4时,h=1/4,有:∫[0, 1] x² dx ≈ (1/4) / 2 * [ f(0) + 2f(1/4) + 2f(2/4) + 2f(3/4) + f(1) ]=(1/4)/2*[0+2*(1/4)²+2*(2/4)²+2*(3/4)²+1]如果我们增大n的取值,将区间分割得更加细致,那么计算得到的近似值将会更加精确。
总结起来,利用复合梯形公式进行积分计算的核心思想是将被积函数分割为若干个小区间,计算每个小区间的积分值,并将它们加总,从而近似得到整个区间上的积分值。
武汉工程大学计算机科学与工程学院计算方法》实验报告日期:年月日实验内容设计分析复化数值积分:将区间[a,b]n 等分,取等距节点x i a ih,iba0,1,2,..., n, hn 由定积分的区间可加性,有nf x dxxix i 1f x dx在每个小区间上利用梯形积分公式有x ii f x dx x i 1 h f x i 1 f x i 2i 1 ih n1T n f a f b 2 f x i一般记 2 i1 称做n+1 点复化梯形积分公式。
数学公式:b h n1f xdx f a f b 2 f x ia 2 i1算法描述:i1Step1:输入a,b 和正整数n;Step2:置h=(b-a)/n;Step3:F=f(a)+f(b);l=0;Step4:对j=1,2,⋯,n 循环执行5;Step5:置x=a+jh; l+=f(x);Step6:置T=h(F+2l)/2 Step7:输出T;程序源代码:#include<iostream>#include<math.h>using namespacestd;<<endl; cout<<" 请输入把 0 到 1 的范围几等分?int m1; cin>>m1;EchelonIntegral(m1); cout<<endl; char answer1;cout<<" 是否要继续求该算法? (y/n)"<<"\t"; cin>>answer1;while(answer1=='y'){cout<<"请输入把 0到 1的范围几等分? "<<"\t"; cin>>m1;EchelonIntegral(m1); //3. 直线求积分; cout<<endl;cout<<" 是否要继续求该算法? (y/n)"<<"\t"; cin>>answer1; }cout<<endl;}cout<<"用梯形积分公式求积分1/(1+pow(sin(x),2))的值测试用例实验总结复化数值积分就是为了减少数值积分的误差,可以把积分区间分成若干小区间,在每个小区间上采取低阶数值积分公式,然后把这些小区间上的数值积分结果加起来作为函数在整个区间上的近似,类似于分段差值。
武汉工程大学计算机科学与工程学院《计算方法》实验报告实验项目复化梯形公式求积分实验类别综合实验实验目的及要求实验目的:学会用复化梯形公式求函数的积分,并应用该算法于实际问题。
实验要求:要求能随意输入被积函数,进行算法设计,打印出误差限例题:求被积函数在0<x<1上的积分))(sin1/(12x+公式:复化梯形公式:设,(i=0,1,…,n-1)nabh-=hxx ii=-+1])(2)()([211∑-=+++=niihafbfafhI误差限:)(12)(//13ininnfhTIfEξ∑-=-=-=成绩评定表类别评分标准分值得分合计上机表现积极出勤、遵守纪律主动完成设计任务30分程序代码比较规范、基本正确功能达到实验要求30分实验报告及时递交、填写规范内容完整、体现收获40分说明:评阅教师:日 期: 年 月 日实 验 内 容设计分析复化数值积分:将区间[a,b]n 等分,取等距节点n ab h n i ih a x i -==+=,,...,2,1,0,由定积分的区间可加性,有()()dxx f dx x f bani x x ii ⎰∑⎰=-=11在每一个小区间上利用梯形积分公式有()()()[]i i x x x f x f hdx x f ii +≈-⎰-112一般记()()()⎥⎦⎤⎢⎣⎡++=∑-=1122n i i n x f b f a f h T 称做n+1点复化梯形积分公式。
数学公式:()()()()⎥⎦⎤⎢⎣⎡++≈∑⎰-=1122n i i bax f b f a f h dx x f 算法描述:Step1:输入a,b 和正整数n ;Step2:置h=(b-a)/n;Step3:F=f(a)+f(b);l=0;Step4:对j=1,2,…,n 循环执行5;Step5:置x=a+jh; l+=f(x);Step6:置T=h(F+2l)/2Step7:输出T;程序源代码:#include<iostream>#include<math.h>using namespace std;double f(double x)//求函数的值;{return1/(1+pow(sin(x),2.0));}void EchelonIntegral(int n)//梯形积分{double y=0;double h=(1.0-0.0)/n;for(int i=0;i<=n;i++){double a=0.0+i*h,b=0.0+(i+1)*h;y=y+h*(f(a)+f(b))/2.0;}cout<<"对应所求的梯形求积分为"<<y<<endl;double En=0.0;double mid=(0.0+1.0)/2.0;double x=mid,p=2*sin(2*x)*sin(2*x)*(1+sin(x)*sin(x))-2*cos(2*x)*pow(1+sin(x)*sin(x),2);double f2d=p/pow(1+sin(x)*sin(x),2);for(i=0;i<n;i++){En=En+pow(h,3)/12.0*f2d;}cout<<"误差为"<<En<<endl;}/*void ParabolicIntegral(int n)//抛物线积分{double y=0;double t=(1.0-0.0)/n;for(int i=0;i<n;i++){double a=0.0+i*t,b=0.0+(i+1)*t;y=y+t*(f(a)+f(b)+4*f((a+b)/2.0))/6.0;}cout<<"对应所求的抛物线求积分为"<<y<<endl;}*/void main(){cout<<"*********************用梯形积分公式求积分1/(1+pow(sin(x),2))的值****************"<<endl;cout<<"请输入把0到1的范围几等分?"<<"\t";int m1;cin>>m1;EchelonIntegral(m1);cout<<endl;char answer1;cout<<"是否要继续求该算法?(y/n)"<<"\t";cin>>answer1;while(answer1=='y'){cout<<"请输入把0到1的范围几等分?"<<"\t";cin>>m1;EchelonIntegral(m1);//3.直线求积分;cout<<endl;cout<<"是否要继续求该算法?(y/n)"<<"\t";cin>>answer1;}cout<<endl;}计算机科学与工程学院测试用例计算机科学与工程学院实验总结复化数值积分就是为了减少数值积分的误差,可以把积分区间分成若干小区间,在每个小区间上采取低阶数值积分公式,然后把这些小区间上的数值积分结果加起来作为函数在整个区间上的近似,类似于分段差值。
用复合梯形求积公式求积实验报告学院:计算机学院班级:信计0902班姓名:***学号:**********1、实验目的1)观察复化梯形公式和复化辛普森公式随区间数n 增加各自误差的减少规律;研究广义积分的数值计算如何将其转化为普通积分,再由已有数值积分方法进行计算;2)利用复化梯形公式和复化辛普森公式计算定积分,编程实现。
2、实验内容分别用复化梯形公式和复化辛普森公式计算定积分,研究随着n 增加各自误差的减少规律dx e I x ⎰+=201,取n=2,4,8,16,精确值为I=4.006994。
3、实验原理原理:将区间[a,b]等分成N 个子区间[x(k),x(k+1)](k=0,1,…..,N -1),h=(b-a )/N ,在每个子区间[x(k),x(k+1)]上使用体形公式,可求得结果。
4、设计思想用上述的的公式来计算dx x x I ⎰=10ln ,然后再将积分区间[0,1]进行n 等分,得到n 个小区间[x i ,x i+1](i=0,1,2,…,n-1),每个小区间的长度为h=(b-a)/n,其中x i =a+i*h 。
按∑⎰-=++==101)]()([2)(n i i i ba x f x f h dx x f T ∑-=++=10)()]()([2n i i x f h b f a f h 进行计算。
5、程序框图6、对应程序#include<iostream.h>#include<math.h>double fun1(double a,double h,int n){int i;double sum=0,m;for(i=1;i<=n-1;i++){m=a+i*h;if(m==0)continue;sum=sum+sqrt(m)*log(m);}return sum;}double get(double a){if(a==0)a=0.00001;return sqrt(a)*log(a);}void main(){double a=0,b=1;int n;while(1){cout<<"请输入要等分的子区间个数(n>0):";cin>>n;double h;double T,f0,f1,f2;h=(b-a)/n;f0=get(a);f1=fun1(a,h,n);f2=get(b);T=h/2*(f0+2*f1+f2);cout<<"步长h为:"<<h<<endl;cout<<"复化梯形公式计算结果T="<<T<<endl;}}7、实验结果8、实验体会通过本次实验我熟悉了用复化梯形公式求数值积分的全过程,并更加熟悉的掌握了复化梯形公式原理。
1 实验6 复合梯形公式 (重点)编写函数,实现复合梯形公式,并求解⎰+1
0411
dx x .
算法:
第一步:输入b a ,,ε,置1←n ,a b h -← 第二步:计算)]()([**5.0b f a f h T +=
第三步:计算])*5.0(*[*5.01
∑-=+++=n j jh h a f h T S 第四步:计算T S e -=,置S T ←,n n 2←,h h *5.0← 第五步:若ε≥e 则转到第三步,否则输出T ,结束。
function u=fhtx(a,b,e)
%单独定义函数,所以这里不出现函数的参数。
end
附分组名单:
组 长
第 1组 冯怡纯 罗 佳 马枨煊 李雄伟 第 2组 葛 韵 屈乾宇 王振海 王小曼 第 3组 魏 玲 姚玉婷 蔡文博 薛 姣 第 4组 吴 萍 郭琦慧 袁 鹏 吴昌磊 第 5组 刘鸿亮 唐 朝 裴成黎 张 宏 第 6组 王晓林 李 维 邱建宇 杨 超 第 7组 陈邯蕾 宋 平 邓春燕 张玲玲 第 8组 季海燕 余 雷 林 瑶 谢腾洲。
合工大计算方法试验报告《计算方法》试验报告班级:学号:姓名:实验一、牛顿下山法1实验目的(1)熟悉非线性方程求根简单迭代法,牛顿迭代及牛顿下山法(2)能编程实现简单迭代法,牛顿迭代及牛顿下山法(3)认识选择迭代格式的重要性(4)对迭代速度建立感性的认识;分析实验结果体会初值对迭代的影响2实验内容(1)用牛顿下山法解方程(初值为0.6)输入:初值,误差限,迭代最大次数,下山最大次数输出:近似根各步下山因子(2)设方程f(x)=x-3x–1=0有三个实根x=1.8793,x=-0.34727,x=-1.53209现采用下面六种不同计算格式,求f(x)=0的根x或xx=;x=;x=;x=;x=;x=x-输入:初值,误差限,迭代最大次数输出:近似根、实际迭代次数3算法基本原理求非线性方程组的解是科学计算常遇到的问题,有很多实际背景.各种算法层出不穷,其中迭代是主流算法。
只有建立有效的迭代格式,迭代数列才可以收敛于所求的根。
因此设计算法之前,对于一般迭代进行收敛性的判断是至关重要的。
牛顿法也叫切线法,是迭代算法中典型方法,只要初值选取适当,在单根附近,牛顿法收敛速度很快,初值对于牛顿迭代至关重要。
当初值选取不当可以采用牛顿下山算法进行纠正。
一般迭代:牛顿公式:牛顿下山公式:图3.1一般迭代算法流程图下山因子下山条件4算法描述一般迭代算法见流程图图3.2牛顿下山算法流程图牛顿下山算法见流程图:5、代码:#include<iostream>#include<fstream>#include<cmath>usingnam espacestd;classsrrt{private:intn;double*a,*xr,*xi;public:srrt(in tnn){n=nn;a=newdouble[n+1];//动态分配内存xr=newdouble[n];xi=newdouble[n];}voidinput();//由文件读入代数方程系数voidsrrt_root();//执行牛顿下山法voidoutput();//输出根到文件并显示~srrt(){delete[]a,xr,xi;}};voidsrrt::input()//由文件读入代数方程系数{inti;charstr1[20];cout<<“\n输入文件名:“;cin>>str1;ifstreamfin(str1);if(!fin){cout<<“\n不能打开这个文件“<<str1<<endl;exit(1);}for(i=n;i>=0;i--)fin>>a[i];//读入代数方程系数fin.close();}voidsrrt::srrt_root()//执行牛顿下山法{intm,i,jt,k,is,it;doublet,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;doubleg ,u,v,pq,g1,u1,v1;m=n;while((m>0)&&(fabs(a[m])+1.0==1.0))m=m-1;if(m<=0){cout<<“\n程序工作失败!“<<endl;return;}for(i=0;i<=m;i++)a[i]=a[i]/a[m];for(i=0;i<=m/2; i++){w=a[i];a[i]=a[m-i];a[m-i]=w;}k=m;is=0;w=1.0;jt=1;while(jt==1){pq=fabs(a[k]);while(pq<1. 0e-12){xr[k-1]=0.0;xi[k-1]=0.0;k=k-1;if(k==1){xr[0]=-a[1]*w/a[0];xi[0]=0.0;return;}pq=fabs(a[k]);}q=log(pq);q=q/(1.0* k);q=exp(q);p=q;w=w*p;for(i=1;i<=k;i++){a[i]=a[i]/q;q=q*p;}x=0.0 001;x1=x;y=0.2;y1=y;dx=1.0;g=1.0e+37;l40:u=a[0];v=0.0;for(i=1;i< =k;i++){p=u*x1;q=v*y1;pq=(u+v)*(x1+y1);u=p-q+a[i];v=pq-p-q;}g1=u*u+v*v;if(g1>=g){if(is!=0){it=1;if(it==0){is=1;dd=sqrt(dx *dx+dy*dy);if(dd>1.0)dd=1.0;dc=6.28/(4.5*k);c=0.0;}while(1==1){c =c+dc;dx=dd*cos(c);dy=dd*sin(c);x1=x+dx;y1=y+dy;if(c<=6.29){it=0 ;break;}dd=dd/1.67;if(dd<=1.0e-07){it=1;break;}c=0.0;}if(it==0)gotol40;}else{it=1;while(it==1){ t=t/1.67;it=0;x1=x-t*dx;y1=y-t*dy;if(k>=50){p=sqrt(x1*x1+y1*y1);q=exp(85.0/k);if(p>=q)it=1;}} if(t>=1.0e-03)gotol40;if(g>1.0e-18){it=0;if(it==0){is=1;dd=sqrt(dx*dx+dy*dy);if(dd>1.0)dd=1.0;dc =6.28/(4.5*k);c=0.0;}while(1==1){c=c+dc;dx=dd*cos(c);dy=dd*sin(c );x1=x+dx;y1=y+dy;if(c<=6.29){it=0;break;}dd=dd/1.67;if(dd<=1.0e -07){it=1;break;}c=0.0;}if(it==0)gotol40;}}if(fabs(y)<=1.0e-06){p=-x;y=0.0;q=0.0;}else{p=-2.0*x;q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for(i=1;i<=k;i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w;xi[k-1]=y*w;k=k-1;if(k==1){xr[0]=-a[1]*w/a[0];xi[0]=0.0;}}else{g=g1;x=x1;y=y1;is=0;if(g<=1.0e-22){if(fabs(y)<=1.0e-06){p=-x;y=0.0;q=0.0;}else{p=-2.0*x;q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for(i=1;i<=k;i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w;xi[k-1]=y*w;k=k-1;if(k==1){xr[0]=-a[1]*w/a[0];xi[0]=0.0;}}else{u1=k*a[0];v1=0.0;for(i=2;i<=k;i++){ p=u1*x;q=v1*y;pq=(u1+v1)*(x+y);u1=p-q+(k-i+1)*a[i-1];v1=pq-p-q;}p=u1*u1+v1*v1;if(p<=1.0e-20){it=0;if(it==0){is=1;dd=sqrt(dx*dx+dy*dy);if(dd>1.0)dd=1.0;dc =6.28/(4.5*k);c=0.0;}while(1==1){c=c+dc;dx=dd*cos(c);dy=dd*sin(c );x1=x+dx;y1=y+dy;if(c<=6.29){it=0;break;}dd=dd/1.67;if(dd<=1.0e -07){it=1;break;}c=0.0;}if(it==0)gotol40;if(fabs(y)<=1.0e-06){p=-x;y=0.0;q=0.0;}else{p=-2.0*x;q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for(i=1;i<=k;i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w;xi[k-1]=y*w;k=k-1;if(k==1){xr[0]=-a[1]*w/a[0];xi[0]=0.0;}}else{dx=(u*u1+v*v1)/p;dy=(u1*v-v1*u)/p;t=1.0+4.0/k;it=1;while(it==1){t=t/1.67;it=0;x1=x-t*dx;y1=y-t*dy;if(k>=50){p=sqrt(x1*x1+y1*y1);q=exp(85.0/k);if(p>=q)it=1;}} if(t>=1.0e-03)gotol40;if(g>1.0e-18){it=0;if(it==0){is=1;dd=sqrt(dx*dx+dy*dy);if(dd>1.0)dd=1.0;dc =6.28/(4.5*k);c=0.0;}while(1==1){c=c+dc;dx=dd*cos(c);dy=dd*sin(c );x1=x+dx;y1=y+dy;if(c<=6.29){it=0;break;}dd=dd/1.67;if(dd<=1.0e -07){it=1;break;}c=0.0;}if(it==0)gotol40;}if(fabs(y)<=1.0e-06){p=-x;y=0.0;q=0.0;}else{p=-2.0*x;q=x*x+y*y;xr[k-1]=x*w;xi[k-1]=-y*w;k=k-1;}for(i=1;i<=k;i++){a[i]=a[i]-a[i-1]*p;a[i+1]=a[i+1]-a[i-1]*q;}xr[k-1]=x*w;xi[k-1]=y*w;k=k-1;if(k==1){xr[0]=-a[1]*w/a[0];xi[0]=0.0;}}}}if(k==1)jt=0;elsejt=1;}}voidsrrt::outp ut()//输出根到文件并显示{intk;charstr2[20];cout<<“\n输出文件名:“;cin>>str2;ofstreamfout(str2);if(!fout){cout<<“\n不能打开这个文件“<<str2<<endl;exit(1);}for(k=0;k<n;k++){fout<<xr[k]<<““<<xi[k ]<<endl;cout<<xr[k]<<“+j“<<x i[k]<<endl;}fout.close();}voidmain ()//主函数{srrtroot(6);root.input();//由文件读入代数方程系数root.srrt_root();//执行牛顿下山法root.output();//输出根到文件并显示}6、输入输出输出结果如下:7、分析体会牛顿下山法作为计算方法课程中重要的知识点,在看书学习时较易大致理解其思路,但上级编写代码时却是有难度的。
数值分析实验 三班级:10信计2班 学号:59 姓名:王志桃 分数一·问题提出:选用复合梯形公式,复合Simpson 公式,计算(1) I =dx x ⎰-4102sin 4 ()5343916.1≈I(2) I = dx x x⎰1sin ()9460831.0,1)0(≈=I f(3) I = dx xe x⎰+1024(4) I = ()dx x x ⎰++10211ln二·实验要求:1.编制数值积分算法的程序2.分别用两种算法计算同一个积分,并比较计算结果3.分别取不同步长()/ a b h -=n ,试比较计算结果(如n = 10, 20等)4.给定精度要求ε,试用变步长算法,确定最佳步长三·实验流程图:复化梯形公式:输入 端点 a , b 正整数 n直接计算TN=h/2*[f(a)+2∑f(x k )+f(b)] k=1,2…,n-1输出 定积分近似值TN复化Simpson 公式输入 端点 a , b 正整数 n输出 定积分近似值SN(1) 置h=(b-a)/(2n)(2) F0=f(a)+f(b) , F1=0 , F2=0(3) 对j=1,2,…,2n-1循环执行步4到步5(4) 置x=a+jh(5) 如果j 是偶数,则F2=F2+f(x),否则F1=F1+f(x)(6) 置SN=h(F0+4F1+2F2)/3(7) 输出SN,停机四·源程序:#include<iostream>#include<math.h>using namespace std;#define n 20//此为步长double f1(double x){double y;y=sqrt(4-sin(x)*sin(x));return y;}double f2(double x){if(x==0)return 1;double y;y=sin(x)/x;return y;}double f3(double x){double y;y=exp(x)/(4+x*x);return y;}double f4(double x){double y;y=log(1+x)/(1+x*x);return y;}int main(){int j;double e=0.000001,h,F0,F1,F2,a,b,x,S;cout<<"利用复化Simpson公式求积分"<<endl;//1a=0;b=0.25*3.141592;h=(b-a)/(2*n);F0=f1(a)+f1(b);F1=F2=0;for(j=1;j<2*n;j++){x=a+j*h;if(j%2==0)F2=F2+f1(x);elseF1=F1+f1(x);}S=((F0+F1*4+F2*2)*h)/3;cout<<"第一个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl;//2a=0;b=1;h=(b-a)/(2*n);F0=f2(a)+f2(b);F1=F2=0;for(j=1;j<2*n;j++){x=a+j*h;if(j%2==0)F2=F2+f2(x);elseF1=F1+f2(x);}S=(F0+F1*4+F2*2)*h/3;cout<<"第二个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl;//3a=0;b=1;h=(b-a)/(2*n);F0=f3(a)+f3(b);F1=F2=0;for(j=1;j<2*n;j++){x=a+j*h;if(j%2==0)F2=F2+f3(x);elseF1=F1+f3(x);}S=(F0+F1*4+F2*2)*h/3;cout<<"第三个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl;//4a=0;b=1;h=(b-a)/(2*n);F0=f4(a)+f4(b);F1=F2=0;for(j=1;j<2*n;j++){x=a+j*h;if(j%2==0)F2=F2+f4(x);elseF1=F1+f4(x);}S=(F0+F1*4+F2*2)*h/3;cout<<"第四个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl<<endl;cout<<"利用复化梯形公式求积分"<<endl;//1a=0;b=0.25*3.141592;h=(b-a)/n;F0=f1(a)+f1(b);F1=0;for(j=1;j<n;j++){x=a+j*h;F1=F1+f1(x);}S=((F0+F1*2)*h)/2;cout<<"第一个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl;//2a=0;b=1;h=(b-a)/n;F0=f2(a)+f2(b);F1=0;for(j=1;j<n;j++){x=a+j*h;F1=F1+f2(x);}S=((F0+F1*2)*h)/2;cout<<"第二个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl;//3a=0;b=1;h=(b-a)/n;F0=f3(a)+f3(b);F1=0;for(j=1;j<n;j++){x=a+j*h;F1=F1+f3(x);}S=((F0+F1*2)*h)/2;cout<<"第三个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl;//4a=0;b=1;h=(b-a)/n;F0=f4(a)+f4(b);F1=0;for(j=1;j<n;j++){x=a+j*h;F1=F1+f4(x);}S=((F0+F1*2)*h)/2;cout<<"第四个积分公式:端点a为"<<a<<"、b为"<<b<<",n为"<<n<<endl<<"结果为"<<S<<endl;return 0;}五.实验结果六.实验心得:通过本次实验,我掌握了求数值积分的各种方法。
合肥工业大学计算方法复化梯形公式实验《计算方法》实验报告学号姓名班级实验项目名称一、实验名称实验二数值积分实验二数值积分二、实验目的:(1)悉复化梯形方法、复化Simpson熟方法、梯形递推算法、龙贝格算法;(2)编程实现复化梯形方法、复化Simpson能方法、梯形递推算法、龙贝格算法;(3)理解并掌握自适应算法和收敛加速算法的基本思想;(4)分析实验结果体会各种方法的精确度,建立计算机求解定积分问题的感性认识三、实验内容及要求(1)设计复化梯形公式求积算法,编制并调试相应的函数子程序(2)设计复化辛浦生求积算法,编制并调试相应的函数子程序( 3)用龙贝格算法计算0sin x dx1x计算机科学与工程学院输入:积分区间,误差限输出:序列 Tn ,Sn,Cn,Rn 及积分结果(参考书本P81 的表 2-5)取n=2,4,8,16,精确解为 0.9460831四、实验原理及算法描述在许多实际问题中,常常需要计算定积分bf (x)dx 的a值。
根据微积分学基本定理,若被积函数f(x) 在区间[a,b] 上连续,只要能找到f(x) 的一个原函数F(x) ,便可利用牛顿 -莱布尼兹公式 b f (x) F (b) F (a) 求得积分值。
a但是在实际使用中,往往遇到如下困难,而不能使用牛顿 -莱布尼兹公式。
(1)找不到用初等函数表示的原函数(2)虽然找到了原函数,但因表达式过于复杂而不便计算(3) f(x) 是由测量或计算得到的表格函数由于以上种种困难,有必要研究积分的数值计算问题。
利用插值多项式P则积分bn ( x ) f ( x )f ( x)dx 转化为a b P n( x )dx,a算机科学与工程学院然易算。
baPn (x)dx称插型求公式。
最的插型求公式是梯形公式和 Simpson 公式,。
当求点提供多,可以分段使用少点的梯形公式和Simpson 公式,并称复化梯形公式、复化 Simpson公式。
如步未知,可以通差限的控制用区逐次分半的策略自取步的方法称自适算法。
梯形推公式出了区分半前后的推关系。
由梯形推公式求得梯形序列,相序列作性合得Simpson 序列 , Simpson 序列作性合得柯特斯序列 , 柯特斯序列作性合的格序列。
若 |R2-R1|< ,出 R2; 否⋯依此推。
如此加工数据的程叫格算法,如下所示:复化梯形公式 T n1n1h[ f ( x0 )2 f ( x k ) f ( x n )]2k1复化 Simpson 公式Sn hn1n1[ f ( x0 ) 2 f ( x k ) f ( x n ) 4 f ( x1 )]6k1k0k2梯形推公式 T2 n T n h n1f ( x k 1) 2 2 k02加平均公式:4T2 n T nS n42 S2 n S n43 C 2n C n4 1421 C n431R n格算法大大加快了差收的速度,由梯形序列O(h2) 提高到格序列的 O(h8)五、程序代码及实验结果1.主程序计算机科学与工程学院int main(){//cout << longbeige(0, 1, 0.0000001);cout<< " 请输入你的区间和误差限度: " << endl;double x, y, z;cin>> x>> y >> z;cout<< " 根据龙贝格算法求出的精确值为: "<< longbeige(x, y, z)<< endl;cout<< "K"<< ""<< "T"<< ""<< "S" << ""<< "C"<< ""<< "R" << endl;cout<< "0"<< ""<<tixing(0, 1, 1)<< endl;cout<< "1"<< ""<<tixing(0, 1, 2)<<""<< xingbusheng(0, 1, 1)<< endl;cout<< "2"<< ""<< tixing(0, 1, 4)<< ""<< xingbusheng(0, 1, 2)<< ""<<(( double (16/15))* xingbusheng(0, 1, 2)) -(double (1/15))*xingbusheng(0, 1, 1)<< endl;cout<< "3"<< ""<<tixing(0,1,8)<<""<<xingbusheng(0,1,4)<<""<< ((double (16 / 15))* xingbusheng(0,1,4)) -( double (1 / 15))*xingbusheng(0,1, 2) <<" "<< longbeige(0,1,0.000000001) <<endl;cout<< "4"<< ""<< tixing(0, 1, 16)<< " "<< xingbusheng(0, 1, 8)<< " "<< (( double (16 / 15))* xingbusheng(0, 1, 4)) - (double (1 / 15))*xingbusheng(0, 1, 2)<< " "<< longbeige(0, 1, 0.00000001)<< endl;;cout<< "the result "<<4<<" is "<< longbeige(0, 1, 0.00000001);return 0;}2.复化梯形公式子程序 :double tixing( double a, double b, int n){double fa = f( a);double fb = f( b);double h = ( b- a) /n;double fxk = 0.0;for ( int k = 1; k <=n - 1; k++){double xk = a + k*h;fxk = f(xk) + fxk;}double res = (h / 2.0)*(fa + 2.0*fxk + fb);return res;}计算机科学与工程学院3.复化辛浦生公式子程序:double xingbusheng(double a, double b, int n){double h = ( b - a) /n;double fa = f(a);double fb = f(b);//double s = fb - fa;double x =a;double fxk12 = 0.0;double fxk = 0.0;/*for (int k = 1; k <=n; k++){x = x + h / (2.0);double fx = f(x);s = s + 4.0*f(x);x = x + h / (2.0);s = s + 2.0*f(x);}s = (h / 6.0)*s;return s;*/for ( int k = 0; k <=n - 1; k++){double xk1 = a + k*h;x = xk1 + h / 2.0;fxk12 = fxk12 + f(x);}for ( int k = 1; k <=n - 1; k++){double xk = a + k*h;fxk = fxk + f(xk);}double s = (h / 6)*(fa + 4 * fxk12 + 2 * fxk + fb);return s;}4.龙贝格公式子程序:double longbeige( double a, double b, double wuchaxian )double h = b - a;double T1 = (h / 2.0)*(f(a) + f( b));int k = 1;double S,x,T2,S2,S1=0,C1=0,C2,R1=0,R2;double Tck[100],Sck[100],Cck[100],Rck[100]; loop1:S = 0;x = a + h / 2.0;while (x< b){S = S + f(x);x = x + h;}T2 = (T1 / 2.0) + (h / 2)*S;// 第一个S2 = T2 + (1 / 3)*(T2 - T1);if(k == 1){Tck[k] = T2;k++;h = h / 2;T1 = T2;S1 = S2;goto loop1;}C2 = S2 + (1 / 15)*(S2 - S1);if(k == 2){C1 = C2;Tck[k]=T2;k++;h = h / 2;T1 = T2;S1 = S2;goto loop1;}R2 = C2 + (1 / 63)*(C2 - C1);if(k == 3){Tck [k]= T2;C1 = C2;k++;h = h / 2;T1 = T2;S1 = S2;goto loop1;}if (fabs(R2 - R1) >=wuchaxian ){R1 = R2;C1 = C2;k++;h = h / 2;T1 = T2;S1 = S2;goto loop1;}else{return R2;}}实验结果:图1计算机科学与工程学院六、实验总结1.梯形公式的收敛速度太慢,所以我们才会选择之后的几种公式加快收敛速度。