龙贝格算法 一、问题分析
1.1龙贝格积分题目
要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y(x)=l sin(tx),则给定雨篷的长度后，求所需要平板材料的长度）。
二、方法原理
2.1龙贝格积分原理
龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。
在变步长的过程中探讨梯形法的计算规律。设将求积区间[a ，b]分为n 个等
分，则一共有n+1个等分点,k x a kh =+,0,1,b a
h k n
-== ，n 。这里用n T 表示复化梯形法求得的积分值，其下标n 表示等分数。
先考察下一个字段[1,k k x x +]，其中点()1
12
1
2
k k k x
x x ++=
+，在该子段上二分前后两个积分值
()()112
k k h
T f x f x +=
+???? ()()21124k k k h T f x f x f x ++?
???
=
++?? ???????
显然有下列关系
2112122
k h
T T f x +??
=+
???
将这一关系式关于k 从0到n-1累加求和，即可导出下列递推公式
12102122n n n k k h T T f x -+=?
?=+ ???
∑
需要强调指出的是，上式中的b a
h n
-=
代表二分前的步长，而 12
12k x
a k h +
?
?=++ ??
?
梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计
算量，自然式人们极为关心的。
根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，既有
211
14
n n T T -≈- 将上式移项整理，知
221
1()3
n n n T T T -≈-
由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计
法。
按上式，积分值2n T 的误差大致等于21()3
n n T T -，如果用这个误差值作为2n T 的一种
补偿，可以期望，所得的
()222141333
n n n n n T T T T T T =+
-=- 应当是更好的结果。
按上式，组合得到的近似值T 直接验证，用梯形二分前后的两个积分值n T 和2n T 按
式组合，结果得到辛普生法的积分值n S 。
24133
n n n S T T =-
再考察辛普生法。其截断误差与4h 成正比。因此，若将步长折半，则误差相
应的减至十六分之一。既有
21
16
n n I S I S -≈- 由此得
21611515
n n I S S ≈
- 不难验证，上式右端的值其实就等于n C ，就是说，用辛普生法二分前后的两个积分值
n S 和2n S ，在按上式再做线性组合，结果得到柯特斯法的积分值n C ，既有
2161
1515
n n n C S S ≈
- 重复同样的手续，依据斯科特法的误差公式可进一步导出龙贝格公式
26416363
n n n R C C ≈-
应当注意龙贝格公式已经不属于牛顿—柯特斯公式的范畴。
在步长二分的过程中运用公式加工三次，就能将粗糙的积分值n T 逐步加工成
精度较高的龙贝格n R ，或者说，将收敛缓慢的梯形值序列n T 加工成熟练迅速的龙贝格值序列n R ，这种加速方法称龙贝格算法。
三、算法设计
3.1龙贝格积分算法
就是求出1T ,再走一遍求出2T ，根据12T T 求出1S ,再走一遍求出4T ，根据2T 4T 求出2S ，
根据1S 2S 求出1C ,再走一遍程序求出8T ，根据4T 8T 得出4S ，根据2S 4S 得出2C ，再根据1C 2C 得出1R ,再走一边程序，
得出16T ，根据8T 16T 得出8S ，根据4S 8S 得出4C ，再由2C 4C 得出2R 。再根据1R 2R 相减的绝对值小于其精度。那其中2R 为求出的值。
四、案例分析
4.1龙贝格积分分析
a—积分下限
b—积分上限
n—区间个数
e—积分值要求达到的精度
s—用以存放除积分区间两端点以外的其他各节点函数值的累加和
p—积分区间两端点函数值之和
h—步长值
T1 、T2分别存放二分区间前后梯形积分值
S1 、S2分别存放二分区间前后辛普生积分值
C1 、C2分别存放二分区间前后斯科特积分值
R1 、R2分别存放二分区间前后龙贝格积分值
五、总结
5.1龙贝格积分总结
通过本次试验，了解了龙贝格算法的计算过程，了解了龙贝格公式的计算收敛过程，用变步长的方法，逐步减小步长，反复积分，逐步得到所求积分值满足精度要求。一步步从梯形法的递推到辛普森到柯特斯法，最后到龙贝格，让精度逐步升高。
附录
龙贝格积分：
#include "stdio.h"
#include "math.h"
float l;
float t;
int main(void){
float f(float);
float a,b,e,h,T1=0,T2=0,S1=0,S2=0,C1=0,C2=0,R1=0,R2=0,k,s,x;
int i=0;
printf("\n****************************************\n");
printf("****************龙贝格算法**************\n");
printf("****************************************\n\n");
printf("请输入积分的下限：");
scanf("%f",&a);
printf("\n请输入积分的上限：");
scanf("%f",&b);
printf("\n请输入允许误差：");
scanf("%f",&e);
printf("请输入L：");
scanf("%f",&l);
printf("请输入T：");
scanf("%f",&t);
k=1;
h=b-a;
T1=h*(f(a)+f(b))/2;
printf("----------------------\n");
printf("k T2 S2 C2 R2\n");
printf("%d %10.7f %10.7f %10.7f %10.7f\n",i,T1,S1,C1,R1); do
{
s=0;
x=a+h/2;
while(x { s+=f(x); x+=h; } T2=T1/2+s*h/2; S2=T2+(T2-T1)/3; if (k==1) { k=k+1; h=h/2; T1=T2; S1=S2; } else if (k==2) { C2=S2+(S2-S1)/15;C1=C2; k=k+1; h=h/2; T1=T2; S1=S2; } else if (k==3) { R2=C2+(C2-C1)/63;C2=S2+(S2-S1)/15;C1=C2;k=k+1;h=h/2;T1=T2;S1=S2; } else { C2=S2+(S2-S1)/15; R2=C2+(C2-C1)/63; if (fabs(R2-R1) { printf("%d %10.7f %10.7f %10.7f %10.7f\n",i+1,T2,S2,C2,R2); break; } else { R1=R2; C1=C2; k=k+1; h=h/2; T1=T2; S1=S2; } } i++; printf("%d %10.7f %10.7f %10.7f %10.7f\n",i,T2,S2,C2,R2); }while (1); getchar(); return 0; } float f(float x) { float y=0; if(x==0.0) return 1; y=(float)l*sin(t*x)/x; return y; } 《数据结构与算法》期末复习题 一、选择题。 1．在数据结构中，从逻辑上可以把数据结构分为 C 。 A．动态结构和静态结构B．紧凑结构和非紧凑结构 C．线性结构和非线性结构D．内部结构和外部结构 2．数据结构在计算机内存中的表示是指 A 。 A．数据的存储结构B．数据结构C．数据的逻辑结构D．数据元素之间的关系 3．在数据结构中，与所使用的计算机无关的是数据的 A 结构。 A．逻辑B．存储C．逻辑和存储D．物理 4．在存储数据时，通常不仅要存储各数据元素的值，而且还要存储 C 。 A．数据的处理方法B．数据元素的类型 C．数据元素之间的关系D．数据的存储方法 5．在决定选取何种存储结构时，一般不考虑 A 。 A．各结点的值如何B．结点个数的多少 C．对数据有哪些运算D．所用的编程语言实现这种结构是否方便。 6．以下说法正确的是 D 。 A．数据项是数据的基本单位 B．数据元素是数据的最小单位 C．数据结构是带结构的数据项的集合 D．一些表面上很不相同的数据可以有相同的逻辑结构 7．算法分析的目的是 C ，算法分析的两个主要方面是 A 。 （1）A．找出数据结构的合理性B．研究算法中的输入和输出的关系C．分析算法的效率以求改进C．分析算法的易读性和文档性 （2）A．空间复杂度和时间复杂度B．正确性和简明性 C．可读性和文档性D．数据复杂性和程序复杂性 8．下面程序段的时间复杂度是O(n2) 。 s =0; for( I =0; i 数值分析课程的大作业，教材《数值分析》李乃成.梅立泉 clear clc format long f=input('请输入原函数f=','s'); a=input('积分下限a='); b=input('积分上限b='); eps1=input('精度eps1='); T(1)=double((b-a)/2*(limit(sym(f),findsym(sym(f)),a)+limit(sym(f),findsym(sym(f)),b))); for k=2:4 sum1=0; for i=1:2^(k-2) sum1=sum1+subs(sym(f),findsym(sym(f)),(a+(2*i-1)*(b-a)/2^(k-1))); end T(k)=1/2*T(k-1)+(b-a)/(2^(k-1))*sum1; end for k=1:3 S(k)=T(k+1)+1/(4-1)*(T(k+1)-T(k)); end for k=1:2 C(k)=S(k+1)+1/(4^2-1)*(S(k+1)-S(k)); end R(1)=C(2)+1/(4^3-1)*(C(2)-C(1)); k=3; while 1 T(1)=T(2); T(2)=T(3); T(3)=T(4); sum2=0; for i=1:2^k sum2=sum2+subs(sym(f),findsym(sym(f)),(a+(2*i-1)*(b-a)/2^(k+1))); end T(4)=1/2*T(4)+(b-a)/2^(k+1)*sum2; S(1)=S(2); S(2)=S(3); S(3)=T(4)+1/(4-1)*(T(4)-T(3)); C(1)=C(2); C(2)=S(3)+1/(4^2-1)*(S(3)-S(2)); R(2)=C(2)+1/(4^3-1)*(C(2)-C(1)); if abs(R(2)-R(1)) 数值分析 实 验 报 告 专业：信息与计算科学 班级： 10***班 学号： 1008060**** 姓名： ****** 实验目的： 用龙贝格积分算法进行积分计算。 算法要求： 龙贝格积分利用外推方法，提高了计算精度，加快了收敛速度。 1--4R R R R 1-j 1-j 1-k 1-j k 1-j k j k ，，，，+= ，k=2，3，… 对每一个k ，j 从2做到k ，一直做到|R R 1-k 1-k k k -，，| 小于给定控制精 度时停止计算。 其中： T R h k 1k =，（复化梯形求积公式），2h 1-k k a -b = 程序代码： #include while(t!=n) { t*=2; ++i; } return i; } float t(int n) { float g,h,q=0; if(1==n) { h = (float)fabs(b-a); q = (f(a)+f(b))*h/2; } else { float x = a; g = 0; h = (float)fabs(b-a)*2/n; x = x+h/2; while(x 一、基本算法 1．交换（两量交换借助第三者） 例1、任意读入两个整数，将二者的值交换后输出。 main() {int a,b,t; scanf("%d%d",&a,&b); printf("%d,%d\n",a,b); t=a; a=b; b=t; printf("%d,%d\n",a,b);} 【解析】程序中加粗部分为算法的核心，如同交换两个杯子里的饮料，必须借助第三个空杯子。 假设输入的值分别为3、7，则第一行输出为3，7；第二行输出为7，3。 其中t为中间变量，起到“空杯子”的作用。 注意：三句赋值语句赋值号左右的各量之间的关系！ 【应用】 例2、任意读入三个整数，然后按从小到大的顺序输出。 main() {int a,b,c,t; scanf("%d%d%d",&a,&b,&c); /*以下两个if语句使得a中存放的数最小*/ if(a>b){ t=a; a=b; b=t; } if(a>c){ t=a; a=c; c=t; } /*以下if语句使得b中存放的数次小*/ if(b>c) { t=b; b=c; c=t; } printf("%d,%d,%d\n",a,b,c);} 2．累加 累加算法的要领是形如“s=s+A”的累加式，此式必须出现在循环中才能被反复执行，从而实现累加功能。“A”通常是有规律变化的表达式，s在进入循环前必须获得合适的初值，通常为0。例1、求1+2+3+……+100的和。 main() {int i,s; s=0; i=1; while(i<=100) {s=s+i; /*累加式*/ i=i+1; /*特殊的累加式*/ } printf("1+2+3+...+100=%d\n",s);} 【解析】程序中加粗部分为累加式的典型形式，赋值号左右都出现的变量称为累加器，其中“i = i + 1”为特殊的累加式，每次累加的值为1，这样的累加器又称为计数器。 数值分析与实验课程设计 班级： 姓名： 学号： 08级应用数学《数值分析与实验(实践)》任务书 一、设计目的 通过《数值分析与实验(实践)》实践环节，掌握本门课程的众多数值解法和原理，并通过编写C语言或matlab程序，掌握各种基本算法在计算机中的具体表达方法，并逐一了解它们的优劣、稳定性以及收敛性。在熟练掌握C语言或matlab语言编程的基础上，编写算法和稳定性均佳、通用性强、可读性好，输入输出方便的程序，以解决实际中的一些科学计算问题。 二、设计教学内容 1、数值方法的稳定性； 2、利用牛顿法和割线法程序求出非线性方程的解，并比较它们之间的优 劣； 3、高斯消去法和列主元高斯消去法求解线性方程组； 雅克比法和高斯－赛德尔迭代法解方程组； 4、利用Lagrange插值多项式求未知点的近似值； 5、利用所给数据进行数据的多项式和可转化成多项式形式的函数拟合； 6、编写复化辛卜生公式和龙贝格算法，通过实际计算体会各种方法的精确 度； 7、利用改进Euler方法和四阶Runge-Kutta方法求解初值问题的微分方 程组； 8、利用幂法求矩阵按模最大的特征值及对应特征向量； （8个中选取1个） 三、设计时间 2011—2012学年第1学期：第16周共计一周 教师签名: 2011年12月12日 前言 数值计算方法是一种利用计算机解决数学问题的数值近似解方法，特别是无法用人工过计算器计算的数学问题。数值计算方法常用于矩阵高次代数方程矩阵特征值与特征向量的数值解法，插值法，线性方程组迭代法，函数逼近，数值积分与微分，常微分方程初值问题数值解等。 作为数学与计算机之间的一条通道，数值计算的应用范围已十分广泛，作为用计算机解决实际问题的纽带，数值算法在求解线性方程组，曲线拟合、数值积分、数值微分，迭代方法、插值法、拟合法、最小二乘法等应用广泛。 数值计算方法是和计算机紧密相连的，现代计算机的出现为大规模的数值计算创造了条件，集中而系统的研究适用于计算机的数值方法是十分必要的。数值计算方法是在数值计算实践和理论分析的基础上发展起来的。 通过数值计算方法与实验将有助于我们理解和掌握数值计算方法基本理论和相关软件的掌握，熟练求解一些数学模和运算，并提高我们的编程能力来解决实际问题。 2008-02-18 18:48 【程序1】 题目：有1、2、3、4个数字，能组成多少个互不相同且无重复数字的三位数？都是多少？ 1.程序分析：可填在百位、十位、个位的数字都是1、2、3、4。组成所有的排列后再去 掉不满足条件的排列。 2.程序源代码： main() { int i,j,k; printf("\n"); for(i=1;i<5;i++) ／*以下为三重循环*/ for(j=1;j<5;j++) for (k=1;k<5;k++) { if (i!=k&&i!=j&&j!=k) /*确保i、j、k三位互不相同*/ printf("%d,%d,%d\n",i,j,k); } } ============================================================== 【程序2】 题目：企业发放的奖金根据利润提成。利润(I)低于或等于10万元时，奖金可提10%；利润高 于10万元，低于20万元时，低于10万元的部分按10%提成，高于10万元的部分，可可提 成7.5%；20万到40万之间时，高于20万元的部分，可提成5%；40万到60万之间时高于 40万元的部分，可提成3%；60万到100万之间时，高于60万元的部分，可提成1.5%，高于 100万元时，超过100万元的部分按1%提成，从键盘输入当月利润I，求应发放奖金总数？ 1.程序分析：请利用数轴来分界，定位。注意定义时需把奖金定义成长整型。 2.程序源代码： main() { long int i; int bonus1,bonus2,bonus4,bonus6,bonus10,bonus; scanf("%ld",&i); bonus1=100000*0.1;bonus2=bonus1+100000*0.75; bonus4=bonus2+200000*0.5; bonus6=bonus4+200000*0.3; bonus10=bonus6+400000*0.15; if(i<=100000) 实验三 龙贝格方法 【实验类型】 验证性 【实验学时】 2学时 【实验内容】 1.理解龙贝格方法的基本思路 2.用龙贝格方法设计算法，编程求解一个数值积分的问题。 【实验前的预备知识】 1．计算机基础知识2．熟悉编程基本思想3．熟悉常见数学函数； 【实验方法或步骤】 龙贝格方法的基本思路龙贝格方法是在积分区间逐次二分的过程中，通过 对梯形之值进行加速处理，从而获得高精度的积分值。 1． 龙贝格方法的算法 步骤1 准备初值()f a 和()f b ，用梯形计算公式计算出积分近似值 ()()12b a T f a f b -=+??? ? 步骤2 按区间逐次分半计算梯形公式的积分近似值令 2i b a h -=，0,1,2,...i =计算12102122n n n i i h T T f x -+=??=+ ??? ∑，2i n = 步骤3 按下面的公式积分梯形公式：()223n n n n T T S T -=+ 辛普生公式：()2215n n n n S S C S -=+ 龙贝格公式：()2263n n n n C C R C -=+ 步骤4 精度控制 当2n n R R ε-<，（ε为精度）时，终止计算，并取2n R 为近似值否则将步长折 半，转步骤2。 [实验程序] #include Romberg龙贝格算法实验报告 2017-08-09 课程实验报告 课程名称： 专业班级： CS1306班学号： U201314967 姓名：段沛云指导教师：报 告日期： 计算机科学与技术学院 目录 1 实验目的 (1) 2 实验原理 (1) 3 算法设计与流程框图 (2) 4 源程序 (4) 5 程序运行 (7) 6 结果分析 (7) 7 实验体会 (7) 1 实验目的 掌握Romberg公式的用法，适用范围及精度，熟悉Romberg算法的流程，并能够设计算法计算积分 31 得到结果并输出。 1x 2 实验原理 2.1 取k=0,h=b-a,求T0= 数)。 2.2 求梯形值T0( b-a ),即按递推公式（4.1）计算T0。 k 2 h [f(a)+f(b)]，令1→k,(k记区间[a，b]的二分次2 2.3 求加速值，按公式(4.12)逐个求出T表的第k行其余各元素Tj(k-j) (j=1,2,….k)。 2.4 若|Tk+1-Tk| n-1 11T2n=[Tn+hn∑f(xi+)] 22i=0 1 Sn=T2n+(T2n-Tn) 31 Cn=S2n+(S2n-Sn) 151 Rn=C2n+(C2n-Cn) 63 3 算法设计与流程框图 算法设计：(先假定所求积分二分最大次数次数为20) 3.1 先求T[k][0] 3.2 再由公式T (k)m 4m(k+1)1)=mTm-1-mTm(k-1(k=1,2,) 求T[i][j] 4-14-1 3.3 在求出的同时比较T[k][k]与T[k-1][k-1]的大小，如果二者之差的绝对 值小于1e-5，就停止求T[k][k]；此时的k就是所求的二分次数，而此时的T[k][k]就是最终的结果 3.4 打印出所有的T[i][j]；程序流程图 #include C1=C; } if(i>=2) R=C+(C-C1)/(pow(4,3)-1); if(abs(R-R1)<0.0000001) /*判断是否满足计算精度，满足则退出迭代，不满足则继续*/ break; R1=R; } printf("积分计算数值为\n"); /*输出计算结果*/ printf("%2.7f\n",T); } double fx1(double x) /*编写一个被积函数fx*/ { double fx; fx=1/(1+x); return(fx); } double fx2(double i,double a,double b) /*编写T的高次迭代函数*/ { //double fx1(double z); double m,l; double Tk,Tkk,z; m=pow(2,i); Tk=0; for(l=1;l<=m;l++) { z=a+(2*l-1)*(b-a)/(2*m); Tk+=fx1(z); } Tkk=(b-a)*Tk/(m*2); return(Tkk); 1.定积分近似计算： /*梯形法*/ double integral(double a,double b,long n) { long i;double s,h,x; h=(b-a)/n; s=h*(f(a)+f(b))/2; x=a; for(i=1;i } 3.素数的判断： /*方法一*/ for (t=1,i=2;i 作业三——龙贝格算法的matlab实现程序流程图： 程序源代码： 文件f.m function fx = f(x) if x == 0 fx = 1; else fx = sin(x) / x; end end 文件longbeige.m clc clear all; format long a=input('请输入你要求得积分的下限:'); b=input('请输入你要求得积分的上限:'); e=input('请输入你要求得积分的结束精度:'); k=input('请输入你要求得积分的最大次数:'); fx=@(x)sin(x)/x; lbg(@f,a,b,k,e) 文件lbg.m function lbg(fx,a,b,k,e) T=zeros(k,k); T(1,1)=(b-a)*(1+fx(b))/2; for i=1:k m=0; for j=1:2^(i-1) m=m+fx(a+(2*j-1)*(b-a)/(2^i)); end T(i+1,1)=0.5*T(i,1)+(b-a)*m/2^i; for n=1:i T(i+1,n+1)=(4^n*T(i+1,n)-T(i,n))/(4^n-1); end if abs(T(i+1,i+1)-T(i,i))<=e & i>=4 break; else ; end end for i=1:k if T(i,1)==0 j=i; break; else ; end end if j==k error('所求次数不够或不可积') else ; end T=T(1:j-1,1:j-1) disp('所求的积分值为：') disp(T(j-1,1)) #include if(i>=j) printf("%10.8f\t",t[i][j]); printf("\n"); } } } void main() { printf("please input a= ,b= ,eps= \n"); scanf("%lf %lf %lf",&a,&b,&eps); romberg(a,b); } 运行结果： please input a= ,b= ,eps= 0.0 1.0 0.00000001 k=0 0.92073549 k=1 0.93979328 0.94614588 k=2 0.94451352 0.94608693 0.94608300 k=3 0.94569086 0.94608331 0.94608307 0.94608307 k=4 0.94598503 0.94608309 0.94608307 0.94608307 0.946083 07 Press any key to continue C语言常用算法集合 1.定积分近似计算： /*梯形法*/ double integral(double a,double b,long n) { long i;double s,h,x; h=(b-a)/n; s=h*(f(a)+f(b))/2; x=a; for(i=1;i if(n==1||n==2) *s=1; else{ fib(n-1,&f1); fib(n-2,&f2); *s=f1+f2; } } 3.素数的判断： /*方法一*/ for (t=1,i=2;i 龙贝格算法 一、问题分析 1.1龙贝格积分题目 要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y(x)=l sin(tx),则给定雨篷的长度后，求所需要平板材料的长度）。 二、方法原理 2.1龙贝格积分原理 龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。 在变步长的过程中探讨梯形法的计算规律。设将求积区间[a ，b]分为n 个等 分，则一共有n+1个等分点,k x a kh =+,0,1,b a h k n -==，n 。这里用n T 表示复化梯形法求得的积分值，其下标n 表示等分数。 先考察下一个字段[1,k k x x +]，其中点()1 12 1 2 k k k x x x ++= +，在该子段上二分前后两个积分值 ()()112 k k h T f x f x += +???? ()()21124k k k h T f x f x f x ++? ???=++?? ??????? 显然有下列关系 2112122 k h T T f x +?? =+ ??? 将这一关系式关于k 从0到n-1累加求和，即可导出下列递推公式 1 210 2122n n n k k h T T f x -+=?? =+ ??? ∑ 需要强调指出的是，上式中的b a h n -= 代表二分前的步长，而 12 12k x a k h + ? ?=++ ?? ? 梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计 算量，自然式人们极为关心的。 根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，既有 211 14 n n T T -≈- 将上式移项整理，知 221 1()3 n n n T T T -≈- 由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计 法。 按上式，积分值2n T 的误差大致等于21()3 n n T T -，如果用这个误差值作为2n T 的一种 补偿，可以期望，所得的 ()222141333 n n n n n T T T T T T =+ -=- 应当是更好的结果。 按上式，组合得到的近似值T 直接验证，用梯形二分前后的两个积分值n T 和2n T 按 式组合，结果得到辛普生法的积分值n S 。 24133 n n n S T T =- 再考察辛普生法。其截断误差与4h 成正比。因此，若将步长折半，则误差相 应的减至十六分之一。既有 21 16 n n I S I S -≈- 由此得 实验八数值积分 信息与计算科学金融崔振威201002034031 一、实验目的： 1、掌握数据积分算法设计及程序实现 二、实验内容： 1、p290-1、p301-2 三、实验要求： 主程序： 复合梯形公式： function [I,step,h2] = CombineTraprl(f,a,b,eps) %f 被积函数 %a,b 积分上下限 %eps 精度 %I 积分结果 %step 积分的子区间数 if(nargin ==3) eps=1.0e-4; end n=1; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; h2=(b-a)/n; function [I,step,h] = IntSimpson(f,a,b,type,eps) %type = 1 辛普森公式 %type = 2 辛普森3/8公式 %type = 3 复合辛普森公式 if(type==3 && nargin==4) eps=1.0e-4; %精度为0.0001 end I=0; switch type case 1, I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+... 4*subs(sym(f),findsym(sym(f)),(a+b)/2)+... subs(sym(f),findsym(sym(f)),b)); step=1; case 2, I=((b-a)/8)*(subs(sym(f),findsym(sym(f)),a)+... 3*subs(sym(f),findsym(sym(f)),(2*a+b)/3)+ ... 3*subs(sym(f),findsym(sym(f)),(a+2*b)/3)+subs(sym(f),findsym(sym(f)),b)); step=1; case 3, n=2; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/6)*(subs(sym(f),findsym(sym(f)),x)+... 4*subs(sym(f),findsym(sym(f)),(x+x1)/2)+... subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; end 60.题目:古典问题:有一对兔子,从出生后第3个月 起每个月都生一对兔子,小兔 子长到第三个月后每个月又生一对兔子,假如兔子都不死,问每个月的兔子总 数 为多少? _________________________________________________________________ _ 程序分析:兔子的规律为数列1,1,2,3,5,8,13,21.... _________________________________________________________________ __ 程序源代码: main() { long f1,f2; int i; f1=f2=1; for(i=1;i<=20;i++) { printf("%12ld %12ld",f1,f2); if(i%2==0) printf("\n");/*控制输出,每行四个*/ f1=f1+f2;/*前两个月加起来赋值给第三个月*/ f2=f1+f2;/*前两个月加起来赋值给第三个月*/ } } 上题还可用一维数组处理,you try! 61.题目:判断101-200之间有多少个素数,并输出所有素数。 _________________________________________________________________ _ 1 程序分析:判断素数的方法:用一个数分别去除2到sqrt(这个数),如果能被 整 除,则表明此数不是素数,反之是素数。 _________________________________________________________________ __ 程序源代码: #include "math.h" main() { int m,i,k,h=0,leap=1; 龙贝格算法应用 目录 一(试验名称 (2) 二(试验目的 (2) 三(方法原理 (2) 四(算法设计 (4) 五(实例分析 (6) 六(总结 (6) 七(附录 (9) 1 一(试验名称 龙贝格算法的应用 二(试验目的 要求学生运用龙贝格算法解决实际问题(塑料雨篷曲线满足函数 y(x)=sin(tx),l则给定雨篷的长度后，求所需要平板材料的长度)。 三(方法原理 (1)实际计算中常采用变步长的计算方案，即在步长逐步减半(即步长二分)的过程中，反复利用复化求积公式进行计算，直到所求得的积分值满足精度要求为止。 1.梯形法的递推化 我们在变步长的过程中探讨梯形法的计算规律。设将求积区间[a,b]分为 ba,h,n等分，则一共有n+1个等分点,,k=0,1...,n。这里用表示用xakh,，Tknn 复化梯形法求得的积分值，其下标n表示等分数。 11x，xx先考虑一下子段[]，其中点=,，在该子段上二分前后xx,，，kkk，1kk，122 的两个积分值 h,，，，Tfxfx ，，，，11kk，，，2 h，，2 ,，，Tfxfxfx，，，，211kk，，，k，,,24，， 显然有下列关系 1h,， TTfx 211，，k，222 将这一关系式关于k从0到n-1累加求和，即可导出下列递推公式 n,11h,， TTfx,，21nn，k，222k,0 ba,h,需要强调指出的是，上式中的代表二分前的步长，而 n 1,,xakh ,，，1,,k，22,, 2 2.龙贝格公式 梯形法的算法简单，单精度低，收敛的速度缓慢。 2 根据梯形法的计算公式(21)，积分值的截断误差大致与成正比，因hTn 1此步长减半后误差将减至，既有 4 1,T12n ,14,Tn 将上式移项整理，知 1,,,TTT (25) 1，，22nnn3 与相当接近，就可以保证计算结由此可见，只要二分前后两个积分值TT2nn果的误差很小。 T2n 1TT,按式(25)，积分值的误差大致等于，如果用这个误差值作为TT，， 2nn2n2n3 多种数值积分求积公式的分析比较 吴春晖 （中国海洋大学海洋环境学院山东青岛 266100） 摘要： 对于运用牛顿-莱布尼茨积分公式不能较好解决的定义在区间[a,b]上的可积函数，原函数并不能简单地用初等函数来表达，故需要构造定积分的近似计算公式。在本文中，主要构建了抛物线求积公式及其复化抛物线公式。在对抛物线类的求积公式进行应用检验后，再运用Gauss求积公式，构建Gauss-Laguerre求积公式，对相同的问题进行运用，并比较截断误差。之后再对求积过程进行优化，在限定误差范围的情况下，利用龙贝格算法，对求积加速收敛。 关键词：抛物线求积复化求积Gauss-Laguerre加速收敛 引言： 对于一些较为复杂的函数，在一定的误差要求下，需要通过构造的方式求给定函数的定积分。基本的替代法主要有梯形面积及抛物线近似代替曲边梯形。并通过划分更小的区间，减少截断误差从而提出了复化梯形及抛物线公式。为了提高运算效率，有加速收敛的Richardson外推法和Romberg求积公式。之后，针对节点数固定情况下，提出了Gauss公式，使其获得最大的精度。本文主要研究的是抛物线求积法与Gauss-Laguerre公式。 目录 第一章抛物线求积公式及应用 (3) 1.1抛物线求积公式的算法 (3) 1.2抛物线求积公式的matlab程序 (3) 1.3复化抛物线求积公式的应用 (4) 第二章Gauss-Laguerre求积公式及应用 (5) 2.1 Gauss-Laguerre的算法 (5) 2.2Gauss-Laguerre公式的matlab程序 (5) 2.3 Gauss-Laguerre求积公式的应用 (6) 第三章龙贝格算法与算法优化 (7) 3.1龙贝格算法及程序 (7) 3.2利用龙贝格算法优化求积 (7) 3.3 龙贝格算法的应用 (8) 第四章数值积分的分析总结 (9) 1.(1) a = 0 b = 1 n = 100 i = 1 h = 1 T1 = 1.3591 i = 2 h = 0.50000 T2 = 0.99734 S1 = 0.87674 i = 3 h = 0.25000 T2 = 0.84726 S2 = 0.79723 C1 = 0.79193 i = 4 T2 = 0.78012 S2 = 0.75774 C2 = 0.75511 R1 = 0.75452 h = 0.062500 T2 = 0.74854 S2 = 0.73801 C2 = 0.73669 R2 = 0.73640 h = 0.031250 T2 = 0.73324 S2 = 0.72815 C2 = 0.72749 R2 = 0.72734 h = 0.015625 T2 = 0.72572 S2 = 0.72321 C2 = 0.72289 R2 = 0.72281 h = 0.0078125 T2 = 0.72199 S2 = 0.72075 C2 = 0.72058 R2 = 0.72055 h = 0.0039062 T2 = 0.72013 S2 = 0.71951 C2 = 0.71943 R2 = 0.71941 h = 0.0019531 T2 = 0.71921 S2 = 0.71890 C2 = 0.71886 R2 = 0.71885 h = 9.7656e-04 T2 = 0.71874 S2 = 0.71859 C2 = 0.71857 R2 = 0.71856 h = 4.8828e-04 T2 = 0.71851 S2 = 0.71844 C2 = 0.71843 R2 = 0.71842 h = 2.4414e-04 T2 = 0.71840 S2 = 0.71836 C2 = 0.71835 R2 = 0.71835 h = 1.2207e-04 T2 = 0.71834 S2 = 0.71832 C2 = 0.71832 R2 = 0.71832 h = 6.1035e-05 T2 = 0.71831 S2 = 0.71830 C2 = 0.71830 R2 = 0.71830 h = 3.0518e-05 T2 = 0.71830 S2 = 0.71829 C2 = 0.71829 R2 = 0.71829 h = 1.5259e-05 T2 = 0.71829 S2 = 0.71829 C2 = 0.71829 R2 = 0.71829 h = 7.6294e-06 T2 = 0.71829 S2 = 0.71828 C2 = 0.71828 R2 = 0.71828 tol = 0.036239 h = 3.8147e-06 T2 = 0.71828 S2 = 0.71828 C2 = 0.71828 R2 = 0.71828 tol = 0.036241 h = 1.9073e-06 T2 = 0.71828 S2 = 0.71828 C2 = 0.71828 R2 = 0.71828 tol = 0.036241 h = 9.5367e-07 T2 = 0.71828 S2 = 0.71828 C2 = 0.71828 R2 = 0.71828 tol = 0.036241 h = 4.7684e-07 T2 = 0.71828 S2 = 0.71828 C2 = 0.71828 R2 = 0.71828 tol = 0.036242 h = 2.3842e-07 T2 = 0.71828 S2 = 0.71828 C2 = 0.71828 R2 = 0.71828 tol = 0.036242 h = 1.1921e-07 T2 = 0.71828 S2 = 0.71828 C2 = 0.71828 R2 = 0.71828 tol = 0.036242数据结构与算法C语言版期末复习题
龙贝格积分-matlab通用程序
数值分析—龙贝格算法
非常全的C语言常用算法
数值分析与实验复化辛卜生公式龙贝格算法
C语言经典算法100例(1---30)
数值分析龙贝格实验报告
Romberg龙贝格算法实验报告.
龙贝格积分法C语言(西安交大)
C语言常用算法集合
龙贝格算法的matlab实现
龙贝格求解积分
最新C语言常用算法集合汇总
龙贝格算法
龙贝格公式和辛普森公式和复合梯形公式
c语言经典算法100例
龙贝格算法应用
多种数值积分的分析比较(Gauss 抛物线 龙贝格)
[2]龙贝格积分[数据]
相关主题
文本预览