执行函数为1、使用方法： Step1：在MATLAB命令窗口输入被积函数 2 1 2 t t e dt 。 输入应为：。 Step2：执行函数。输入形式为mymulNewtonCotes(ft,a,b,m,n)； 其中ft—被积函数，此体重，已经在第一步赋值； a—积分下限，本题中为0； b—积分上限，本题中为1； m—将区间[a,b]等分的子区间数量，本题可选为10； n—采用的Newton-Cotes公式的阶数，必须满足n<8，否则积分没法保证稳定性。 当n=1时，即为复化梯形公式；n=2时，即为复化复化辛普森公式。 所以，分别输入mymulNewtonCotes(ft,0,1,10,1)和mymulNewtonCotes(ft,0,1,10,2)就可以得到两种方法的积分计算结果。 2、计算结果 而根据积分运算，可得： 说明复化梯形和复化辛普森公式计算出的结果基本一致，与实际结果相符。 3、程序代码 functionyy=mymulNewtonCotes(ft,a,b,m,n) %复化Newton-Cotes数值积分公式，即在每个子区间上使用Newton-Cotes公式，然后求和, %参考的输入形式为mymulNewtonCotes(ft,0,1,10,2) %参数说明: %ft——被积函数，此题中ft=@(t)t.*exp(t^2/2) %a——积分下限 %b——积分上限 %m——将区间[a,b]等分的子区间数量 %n——采用的Newton-Cotes公式的阶数，必须满足n<8，否则积分没法保证稳定性 %(1)n=1时为复化梯形公式
%(2)n=2时为复化辛普森公式 xx=linspace(a,b,m+1); forl=1:m s(l)=myNewtonCotes(ft,xx(l),xx(l+1),n); end yy=sum(s); function[y,Ck,Ak]=myNewtonCotes(ft,a,b,n) %牛顿-科特斯数值积分公式 %Ck——科特斯系数 %Ak——求积系数 %y——牛顿-科特斯数值积分结果 xk=linspace(a,b,n+1); forj=1:n+1 ff(j)=ft(xk(j)); end %计算科特斯系数 fori=1:n+1 k=i-1; Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl(@(t)intfun(t,n,k),0,n); end %计算求积系数 Ak=(b-a)*Ck; %求和算积分 y=Ak*ff'; functionf=intfun(t,n,k) %科特斯系数中的积分表达式 f=1; fori=[0:k-1,k+1:n] f=f.*(t-i); end
数值分析 实 验 报 告 专业：信息与计算科学 班级： 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 数值计算课程PROJECT 牛顿-柯特斯公式(Newton-Cotes Methods) 例题介绍： 请在MATLAB环境下编写牛顿-柯特斯公式，以求得π的近似值。 Instruction: Using the Newton-Cotes method to solve all the following problems. Solving the problems with MATLAB, printing out your MATLAB code, figures, as well as necessary problem-solving procedures on paper, and submit before final exam. (No late submission is accepted) Students must solve all these questions correctly to get 5 point extra credits, no partial credit is given. Plagiarizing from other’s work will be treated as 0. Problems I: The value of π can be calculated from the integral 1 2121dx x π-=+?(a)Evaluate the integral by using rectangle method, using 60 subintervals. (b)Evaluate the integral by using midpoint method, using 60 subintervals. (c)Evaluate the integral by using trapezoidal method, using 60 subintervals. (d)Evaluate the integral by using Simpson’s 1/3 method, using 60 subintervals. (e)Evaluate the integral by using Simpson’s 3/8 method, using 60 subintervals. (f)Compare the results and discuss the error from each method. 作业三——龙贝格算法的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)) 摘要 在工程实验及研究中，实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系.可以说，曲线拟合模型与我们的生活生产密切相关. 本课题着重介绍曲线拟合模型及其应用，其中包括它的基本思想、模型的建立、以及具体应用.为了更好的了解曲线拟合模型，可以将它分为线性与非线性模型，在模型建立的基础上我们可以用最小二乘法来解决一些我们日常所应用的问题. 关键词曲线拟合；线性与非线性模型；最小二乘发 目录 引言 (1) 第一章曲线拟合 (2) §1.1 基本思想及基本概念 (2) §1.1.1 方法思想 (2) §1.1.2几个基本概念 (2) §1.2辛普森算法基本定义及其应用 (4) §1.2.1辛普森求积公式的定义 (4) §1.2.2辛普森求积公式的几何意义 (5) §1.2.3辛普森求积公式的代数精度及其余项 (5) §1.2.4辛普森公式的应用 (6) 第二章辛普森求积公式的拓展及其应用 (7) §2.1 复化辛普森求积公式 (7) §2.1.1问题的提出 (7) §2.1.2复化辛普森公式及其分析 (7) §2.1.3复化辛普森公式计算流程图 (8) §2.1.4复化辛普森公式的应用 (9) §2.2 变步长辛普森求积公式 (10) §2.2.1变步长辛普森求积公式的导出过程 (10) §2.2.2变步长辛普森求积公式的加速过程 (12) §2.2.3变步长辛普森求积公式的算法流程图 (13) §2.2.4变步长辛普森公式算法程序代码 (14) §2.2.5变步长辛普森求积公式的应用 (14) §2.2.6小结 (14) §2.2.7数值求积公式在实际工程中的应用 (14) 参考文献 (16) 附录A (17) 实验八数值积分 信息与计算科学金融崔振威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 数值分析与实验课程设计 班级: 姓名: 学号: 08级应用数学《数值分析与实验（实践）》任务书 一、设计目的 通过《数值分析与实验（实践）》实践环节，掌握本门课程的众多数值解法和原理，并通过编写C语言或matlab程序，掌握各种基本算法在计算机中的具体表达方法，并逐一了解它们的优劣、稳定性以及收敛性。在熟练掌握C语言或matlab语言编程的基础上，编写算法和稳定性均佳、通用性强、可读性好，输入输出方便的程序，以解决实际中的一些科学计算问题。 二、设计教学内容 1、数值方法的稳定性； 2、禾U用牛顿法和割线法程序求出非线性方程的解，并比较它们之间的优 劣； 3、高斯消去法和列主元高斯消去法求解线性方程组； 雅克比法和高斯-赛德尔迭代法解方程组； 4、利用Lagrange插值多项式求未知点的近似值； 5、利用所给数据进行数据的多项式和可转化成多项式形式的函数拟合； 6、编写复化辛卜生公式和龙贝格算法，通过实际计算体会各种方法的精确 \ 度； 7、利用改进Euler方法和四阶Runge-Kutta方法求解初值问题的微分方程组； &利用幕法求矩阵按模最大的特征值及对应特征向量； \ （8个中选取1个） 二、设计时间 2011 —2012学年第1学期：第16周共计一周 教师签名: 2011年12月12日 、八 刖 数值计算方法是一种利用计算机解决数学 .言 问题的数值近似解方法， 特别是无法用人工过计算器计算的数学问题。数值计算方法常用于矩阵高次代数方程矩阵特征值与特征向量的数值解法，插值法，线性方程组迭代法，函数逼近，数值积分与微分，常微分方程初值问题数值解等。 作为数学与计算机之间的一条通道，数值计算的应用范围已十分广泛，作为用计算机解决实际问题的纽带，数值算法在求解线性方程组，曲线拟合、数值积分、数值微分，迭代方法、插值法、拟合法、最小二乘法等应用广泛。 数值计算方法是和计算机紧密相连的，现代计算机的出现为大规模的数值计 算创造了条件，集中而系统的研究适用于计算机的数值方法是十分必要的。数值计算方法是在数值计算实践和理论分析的基础上发展起来的。 通过数值计算方法与实验将有助于我们理解和掌握数值计算方法基本理论和相关软件的掌握，熟练求解一些数学模和运算，并提高我们的编程能力来解决实际问题。 龙贝格算法 一、问题分析 1、1龙贝格积分题目 要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y（x)=l sin （tx），则给定雨篷得长度后，求所需要平板材料得长度）. 二、方法原理 2、１龙贝格积分原理 龙贝格算法就是由递推算法得来得。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式. 在变步长得过程中探讨梯形法得计算规律.设将求积区间［ａ，b]分为ｎ个等分，则一共有n+1个等分点，n.这里用表示复化梯形法求得得积分值，其下标n 表示等分数。 先考察下一个字段[]，其中点,在该子段上二分前后两个积分值 显然有下列关系 将这一关系式关于k从０到n－１累加求与，即可导出下列递推公式 需要强调指出得就是,上式中得代表二分前得步长，而 梯形法得算法简单，但精度低,收敛速度缓慢,如何提高收敛速度以节省计算量，自然式人们极为关心得. 根据梯形法得误差公式，积分值得截断误差大致与成正比，因此步长减半后误差将减至四分之一,既有 将上式移项整理,知 由此可见，只要二分前后两个积分值与相当接近，就可以保证计算保证结果计算结果得误差很小，这种直接用计算结果来估计误差得方法称作误差得事后估计法。 ?按上式,积分值得误差大致等于，如果用这个误差值作为得一种补偿，可以期望,所得得 应当就是更好得结果。 ?按上式，组合得到得近似值直接验证，用梯形二分前后得两个积分值与按式组合，结果得到辛普生法得积分值。 再考察辛普生法。其截断误差与成正比.因此,若将步长折半,则误差相应得减至十六分之一。既有 由此得 不难验证,上式右端得值其实就等于,就就是说,用辛普生法二分前后得两个积分值与，在按上式再做线性组合,结果得到柯特斯法得积分值，既有 重复同样得手续，依据斯科特法得误差公式可进一步导出龙贝格公式 应当注意龙贝格公式已经不属于牛顿—柯特斯公式得范畴. 在步长二分得过程中运用公式加工三次，就能将粗糙得积分值逐步加工成精度较高得龙贝格,或者说,将收敛缓慢得梯形值序列加工成熟练迅速得龙贝格值序列，这种加速方法称龙贝格算法。 三、算法设计 3、1龙贝格积分算法 就就是求出,再走一遍求出，根据求出,再走一遍求出，根据求出，根据求出，再走一遍程序求出，根据得出，根据得出,再根据得出,再走一边程序，得出，根据得出,根据得出，再由得出。再根据相减得绝对值小于其精度。那其中为求出得值. 执行函数为mymulNewtonCotes.m 1、使用方法： Step1：在MATLAB 命令窗口输入被积函数212 0t t e dt ?。 输入应为：ft=@(t)t.*exp(t^2/2)。 Step2：执行函数。输入形式为mymulNewtonCotes(ft,a,b,m,n)； 其中ft —被积函数，此体重ft=@(t)t.*exp(t^2/2)，已经在第一步赋值； a —积分下限，本题中为0； b —积分上限，本题中为1； m —将区间[a,b]等分的子区间数量，本题可选为10； n —采用的Newton-Cotes 公式的阶数，必须满足n<8，否则积分没法保 证稳定性。 当n=1时，即为复化梯形公式；n=2时，即为复化复化辛普森公式。 所以，分别输入mymulNewtonCotes(ft,0,1,10,1)和 mymulNewtonCotes(ft,0,1,10,2)就可以得到两种方法的积分计算结果。 2、计算结果 而根据积分运算，可得： 221112 110222 220000() 1.648710.64872t t x x t t e dt e d e dx e e e ====-=-=??? 说明复化梯形和复化辛普森公式计算出的结果基本一致，与实际结果相符。 3、程序代码 function yy = mymulNewtonCotes(ft,a,b,m,n) % 复化Newton-Cotes数值积分公式，即在每个子区间上使用Newton-Cotes公式，然后求和, % 参考的输入形式为mymulNewtonCotes(ft,0,1,10,2) % 参数说明: % ft——被积函数，此题中ft=@(t)t.*exp(t^2/2) % a——积分下限 % b——积分上限 % m——将区间[a,b]等分的子区间数量 % n——采用的Newton-Cotes公式的阶数，必须满足n<8，否则积分没法保证稳定性 % (1)n=1时为复化梯形公式 % (2)n=2时为复化辛普森公式 xx = linspace(a,b,m+1); for l = 1:m s(l) = myNewtonCotes(ft,xx(l),xx(l+1),n); end yy = sum(s); function [y,Ck,Ak] = myNewtonCotes(ft,a,b,n) % 牛顿-科特斯数值积分公式 % Ck——科特斯系数 % Ak——求积系数 % y——牛顿-科特斯数值积分结果 xk = linspace(a,b,n+1); for j = 1:n+1 ff(j) = ft(xk(j)); end % 计算科特斯系数 for i=1:n+1 一．介绍 龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。作为一种外推算法, 它在不增加 开始 计算步长h 计算初值 f(a)、f(b)、 R(1,1) R矩阵迭代计算 No 误差达到精度要 求 Yes 输出R(j+1,j+1) 结束三．源码 1. f=inline('1/(x+1)');%输入函数 a=0;b=1;%取值边界 eps=10^(-7); h=b-a; R(1,1)=h*(f(a)+f(b))/2; j=0; err=1; m=1; while err>eps j=j+1; h=h/2; S=0; for i=1:m x=a+h*(2*i-1); S=S+f(x); end m=2*m; R(j+1,1)=R(j,1)/2+h*S; for i=1:j R(j+1,i+1)=R(j+1,i)+(R(j+1,i)-R(j,i))/(4^i-1); end err=abs(R(j+1,j)-R(j+1,j+1)); end ans=vpa(R(j+1,j+1),7) 2. f=inline('log(x+1)/(x^2+1)');%输入函数 a=0;b=1;%取值边界 eps=10^(-7); h=b-a; R(1,1)=h*(f(a)+f(b))/2; j=0; err=1; m=1; while err>eps j=j+1; h=h/2; S=0; for i=1:m x=a+h*(2*i-1); S=S+f(x); end m=2*m; R(j+1,1)=R(j,1)/2+h*S; for i=1:j R(j+1,i+1)=R(j+1,i)+(R(j+1,i)-R(j,i))/(4^i-1); end err=abs(R(j+1,j)-R(j+1,j+1)); end ans=vpa(R(j+1,j+1),7) 3. f=inline('log(x+1)/x');%输入函数 a=0;b=1;%取值边界 eps=10^(-7); h=b-a; R(1,1)=h*(1+log(2))/2; j=0; err=1; m=1; 龙贝格算法 一、问题分析 1.1龙贝格积分题目 要求学生运用龙贝格算法解决实际问题(塑料雨篷曲线满足函数y(x)=lsin(tx),则给定雨篷的长度后，求所需要平板材料的长度)。 二、方法原理 2.1龙贝格积分原理 龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。 在变步长的过程中探讨梯形法的计算规律。设将求积区间[a，b]分为n个等 ba,分，则一共有n+1个等分点，n。这里用表示复化xakh,，,hk,,,0,1,Tnkn 梯形法求得的积分值，其下标n表示等分数。 1xxx,，先考察下一个字段[]，其中点，在该子段上二分前后xx,，，11kk，kk，1k，22 两个积分值 h ,，，，Tfxfx，，，，11kk，，，2 ，，,,hTfxfxfx,，，，，，，,,,,211kk，k，4,,2,,，， 显然有下列关系 ,,1h TTfx,，,,211k，222,, 将这一关系式关于k从0到n-1累加求和，即可导出下列递推公式 n,1,,1h TTfx,，,,,21nnk，22k,02,, ba,需要强调指出的是，上式中的代表二分前的步长，而 h,n 1,,xakh,，， 1,,k，2,,2 梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计算量，自然式人们极为关心的。 2根据梯形法的误差公式，积分值的截断误差大致与成正比，因此步长减Thn 半后误差将减至四分之一，既有 1,T12n, 14,Tn将上式移项整理，知 1 ,,,TTT1()22nnn3由此可见，只要二分前后两个积分值和相当接近，就可以保证计算保证结TTn2n果计算结果的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计T2n 法。 1 按上式，积分值的误差大致等于，如果用这个误差值作为的一种 TTT,T()2n2nn2n3补偿，可以期望，所得的 141 TTTTTT,，,,,，，222nnnnn333应当是更好的结果。 按上式，组合得到的近似值直接验证，用梯形二分前后的两个积分值和按TTTn2n 式组合，结果得到辛普生法的积分值。 Sn 41 STT,,nnn233 4 再考察辛普生法。其截断误差与成正比。因此，若将步长折半，则误差相h 应的减至十六分之一。既有 IS,12n, IS,16n由此得 161 ISS,,2nn1515 不难验证，上式右端的值其实就等于，就是说，用辛普生法二分前后的两个积分值Cn 和，在按上式再做线性组合，结果得到柯特斯法的积分值，既有 SCS2nnn 实验报告 课程名称：《数值计算方法》 实验名称：数值积分 实验类型：验证性■综合性□设计性□ 实验室名称：数学实验室 班级学号： 09072 学生姓名： 任课教师（教师签名）： 成绩： 实验日期：2012年3月29日 一、实验目的及题目 实验目的：掌握利用复化辛普森公式和龙贝格方法计算积分，掌握复化辛普森公式以及龙贝格方法的原理，熟悉matlab的操作。 题目：利用复化辛普森公式和龙贝格方法计算下列积分： 1、dx e x ?-5 .002 2、dx x x ?20 2 sin )2sin(cos π 二、实验原理、程序框图、程序代码等 实验原理： 1、复化求积公式 由于高阶插值的不稳定性，为了提高计算积分的精度，可把积分区间分为若干个小区间， 将()I f 写成这些小区间上的积分之和，然后对每一个小区间上的积分应用到辛普森公式，或柯特斯公式，并把每个小区间上的结果累加，所得到的求积公式就称为复化求积公式。 将求积区间作n 等分，并记,(0),i b a h x a ih i n n -= =+≤≤于是 11 ()()k k n x x k I f f x dx +-== ∑? 并记 1()()k k x k x I f f x dx += ? 1.1复化辛普森公式 记112 1()2 k k k x x x ++ = +，对每一个()k I f 应用辛普森公式，得到复化辛普森公式 11 01102()()2()()4()6n n n k n k k k h S f f x f x f x f x --+==? ?=+++???? ∑∑ 其截断误差为 14 (4) ()()() ()1802n n k k h h I f S f f η-=-=- ∑ 设[]4(),f x C a b ∈，则 (3) (3) 4 1()()[()()]()180 2 n h I f S f f a f b -=- - 1.2、复化柯斯特公式 高斯—勒让德积分公式 摘要： 高斯—勒让德积分公式可以用较少节点数得到高精度的计算结果,是现在现实生活中经常运用到的数值积分法。然而，当积分区间较大时，积分精度并不理想。 T he adva ntage of Gauss-Legendre integral formula is tend to get high-precision calculational result by using fewer Gauss-points, real life is now often applied numerical integration method. But the precision is not good when the length of integral interval is longer. 关键字： 积分计算，积分公式，高斯—勒让德积分公式，MATLAB Keyword: Integral Calculation , Integral formula ,Gauss-Legendre integral formula, Matlab 引言： 众所周知，微积分的两大部分是微分与积分。微分实际上是求一函数的导数，而积分是已知一函数的导数，求这一函数。所以，微分与积分互为逆运算。 实际上，积分还可以分为两部分。第一种，是单纯的积分，也就是已知导数求原函数，称为不定积分。 相对而言，另一种就是定积分了，之所以称其为定积分，是因为它积分后得出的值是确定的，是一个数，而不是一个函数。 计算定积分的方法很多，而高斯—勒让德公式就是其中之一。 高斯积分法是精度最高的插值型数值积分，具有2n+1阶精度，并且高斯积分总是稳定。而高斯求积系数，可以由Lagrange多项式插值系数进行积分得到。 高斯—勒让德求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定让函数f(x)以此取i=0,1,2....n次多项式使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是2n-1，而且是最高的。通常运用的是（-1，1）的积分节点和积分系数，其他积分域是通过变换 龙贝格积分实验报告 Prepared on 24 November 2020 二、R o m b e r g 积分法 1.变步长Romberg 积分法的原理 复化求积方法对于提高精度是行之有效的方法，但复化公式的一个主要缺点在于要事先估计出部长。若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。而用复化公式的截断误差来估计步长，其结果是步长往往过小，而且''()f x 和(4)()f x 在区间[,]a b 上的上界M 的估计是较为困难的。在实际计算中通常采用变步长的方法，即把步长逐次分半（也就是把步长二等分），直到达到某种精度为止，这种方法就是Romberg 积分法的思想。 在步长的逐步分半过程中，要解决两个问题： 1. 在计算出N T 后，如何计算2N T ,即导出2N T 和N T 之间的递推公式； 2. 在计算出N T 后，如何估计其误差，即算法的终止的准则是什么。 首先推导梯形值的递推公式，在计算N T 时，需要计算1N +个点处的函数值在计算出N T 后，在计算2N T 时，需将每个子区间再做二等分，共新增N 个节点。为了避免重复计算，计算2N T 时，将已计算的1N +个点的数值保留下来，只计算新增N 个节点处的值。为此，把2N T 表示成两部分之和，即 由此得到梯形值递推公式 因此 由复化梯形公式的截断误差有 若''()f x 变化不大时，即''''12()()f f ηη≈，则有 式（2）表明，用2N T 作为定积分I 的近似值，其误差大致为21 ()3 N N T T -，因此其 终止条件为 其中ε是预先给定的精度。 积分公式 将上述方法不断推广下去，可以得到一个求积分的序列，而且这个序列很快收敛到所求的定积分。记 (0)N N T T =，将区间N 等分的梯形值。(1)N N T S =，将区间N 等分的Simpson (2)N N T C =，将区间N 等分的Cotes 。(3)N N T R =，将区间N 等分的Romberg 。 由其可构造一个序列(){}k N T ，次序列称为Romberg 序列，并满足如下递推关系： 以上递推公式就是Romberg 积分递推公式。 积分程序 1. 置1N =，精度要求ε，1h b a =-； 2. 计算(0)1[()()]2 b a T f a f b -=+； 3. 置22 N N h h =，并计算(0) (0)211((21))222N N N k b a b a T T f a k N N =--=++-∑； 4. 置,2,1;M N N N K === 5. 计算(1)(1) 2441 k k k k M M M k T T T ---=-； 6. 若 1M =，则转（7）；否则置2 M M = ，1k k =+转（5）； 7. 若()(1)11k k T T ε--≤，则停止计算（输出()1k T ），否则转（3）。 积分法的应用 function [T,n] = romb(f,a,b,eps) double R ; if nargin<4,eps=1e-8; end 数值分析实验报告 实验3 数值积分 实验目的 通过本实验理解数值积分与微分的基本原理。掌握数值积分中常见的复合求积公式的编程实现。 掌握龙贝格算法的基本思路和迭代步骤； 培养编程与上机调试能力。 算法描述 3.2.1 龙贝格算法基本思路 （1）将区间[]b a ,划分为n 等分，分点：()n x x 0；根据梯形公式 ()()()?? ????+-=∑-=1122n k k n b f x f a f h T ，求出n T ，再根据n T 和n T 2之间的递推公式()∑-=++=10 21222n k k n n x f h T T 求出n T 2； （2）设m 为加速次数，k 为划分区间次数，则由加速公式：()()()k m m k m m m k m T T T 111 41144-----=（ ,2,1=k ）求出第k 次划分，第m 次加速次数的梯形值()k m T ，这样不断地循环，直到求出在满足精度条件下的某个()k m T 作为积分值为止。 3.2.2 龙贝格算法计算步骤 1.输入MATLAB 程序 function[t]=romberg(f,a,b,e) t=zeros(15,4); t(1,1)=(b-a)/2*(f(a)+f(b)); for k=2:4 sum=0; for i=1:2^(k-2) sum=sum+f(a+(2*i-1)*(b-a)/2^(k-1)); end t(k,1)=*t(k-1,1)+(b-a)/2^(k-1)*sum; for i=2:k t(k,i)=(4^(i-1)*t(k,i-1)-t(k-1,i-1))/(4^(i-1)-1); end 二、Romberg 积分法 1.变步长Romberg 积分法的原理 复化求积方法对于提高精度是行之有效的方法，但复化公式的一个主要缺点在于要事先估计出部长。若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。而用复化公式的截断误差来估计步长，其结果是步长往往过小，而且''()f x 和(4)()f x 在区间[,]a b 上的上界M 的估计是较为困难的。在实际计算中通常采用变步长的方法，即把步长逐次分半（也就是把步长二等分），直到达到某种精度为止，这种方法就是Romberg 积分法的思想。 在步长的逐步分半过程中，要解决两个问题： 1. 在计算出N T 后，如何计算2N T ,即导出2N T 和N T 之间的递推公式； 2. 在计算出N T 后，如何估计其误差，即算法的终止的准则是什么。 首先推导梯形值的递推公式，在计算N T 时，需要计算1N +个点处的函数值在计算出N T 后，在计算2N T 时，需将每个子区间再做二等分，共新增N 个节点。为了避免重复计算，计算2N T 时，将已计算的1N +个点的数值保留下来，只计算新增N 个节点处的值。为此，把2N T 表示成两部分之和，即 21 2221 122211 122211 1 [()()2()]21 [()()2()2((21))]21[()()2()]((21))]22N N N N k N N N N N k k N N N N N N k k T h f a f b f a kh h f a f b f a kh f a k h h f a f b f a kh h f a k h -=-==-===+++=+++++-=+++++-∑∑∑∑∑ 由此得到梯形值递推公式 2221 1 ((21))12N N N N N k T T h f a k h ==++-∑ 因此 1 11212122,[()()],2 11 ,() 22h h b a T f a f b h h T T h f a h =-= +==++ // 龙贝格积分公式（2）.cpp : 定义控制台应用程序的入口点。 // #include"stdafx.h" #include数值计算课程PROJECT- 牛顿-柯特斯公式(Newton-Cotes Methods)
龙贝格算法的matlab实现
辛普森求积公式
龙贝格公式和辛普森公式和复合梯形公式
数值分析与实验复化辛卜生公式龙贝格算法
龙贝格算法
matlab实现复化Newton-Cotes公式求积分的程序应用和代码
龙贝格积分
龙贝格算法
数值分析 复化求积 龙贝格 实验报告
数值分析 高斯—勒让德积分公式
龙贝格积分实验报告
龙贝格算法
龙贝格积分实验报告
龙贝格积分公式2
相关主题
文本预览