积分方程的数值计算技巧
- 格式:doc
- 大小:147.19 KB
- 文档页数:8
化学反应速率方程式类积分的数值求解化学反应速率方程式是描述反应速率与物质浓度之间关系的方程。
对于复杂的反应机制,速率方程往往是非线性的,并且很难解析地求得其解析解。
因此,数值求解成为研究化学反应动力学的重要工具之一。
本文将介绍化学反应速率方程式类积分的数值求解方法。
化学反应速率方程式描述了反应速率和反应物浓度之间的关系。
在一般情况下,它可以表示为以下形式:\[ \frac{{dC}}{{dt}} = f(C) \]其中,\(C\) 是反应物的浓度,\(t\) 是时间,\(f(C)\) 是关于浓度的函数。
求解该方程,即求解从初始浓度(\(C_0\))到某一特定时间(\(t\))的浓度变化。
由于反应物浓度的变化是与时间相关的,因此我们需要使用数值方法求解该方程。
常见的化学反应速率方程式类积分的数值求解方法包括 Euler方法、Runge-Kutta 方法和 Adams 方法等。
下面将分别介绍这些方法的基本原理及其应用。
Euler 方法是最简单的数值求解方法之一,它基于离散化的思想,通过将时间区间分割为若干个小的时间步长来近似求解。
对于反应速率方程式来说,Euler 方法的基本形式为:\[ C_{n+1} = C_n + \Delta t \cdot f(C_n) \]其中,\(C_{n+1}\) 是下一个时间步长的浓度,\(C_n\) 是当前时间步长的浓度,\(\Delta t\) 是时间步长,\(f(C_n)\) 是在当前浓度下的反应速率。
Euler 方法的优点是简单易懂,但由于其线性插值的近似方法,其稳定性和精度相对较低。
Runge-Kutta 方法是一类更加精确和稳定的数值求解方法,其基本思想是通过在每个时间步长内进行多次函数评估来提高数值解的精度。
其中最常用的是四阶 Runge-Kutta 方法,其形式如下:\[ \begin{align*} k_1 &= \Delta t \cdot f(C_n) \\ k_2 &= \Delta t\cdot f(C_n + \frac{k_1}{2}) \\ k_3 &= \Delta t \cdot f(C_n +\frac{k_2}{2}) \\ k_4 &= \Delta t \cdot f(C_n + k_3) \\ C_{n+1} &=C_n + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} \end{align*} \]其中,\(k_i\) 是函数在当前浓度和补偿项处的评估。
牛顿迭代法数值积分牛顿迭代法(Newton's method)是一种用于求解方程的迭代数值计算方法,通过不断逼近方程的根来获得精确的解。
其基本思想是利用函数在某点的切线来逼近方程的根,然后通过不断迭代计算来逼近真实的根。
具体而言,假设要求解方程f(x) = 0,首先选择一个初始近似解x_0,然后通过切线的斜率来确定下一个近似解x_1。
切线的斜率可以通过函数的导数f'(x) 来计算,即:k = f'(x_0)。
然后,利用直线的斜截式公式y = k(x - x_0) + f(x_0),将其与x 轴相交得到新的近似解x_1,即使得f(x_1) = 0 的解。
迭代过程如下:1. 选择初始近似解x_0。
2. 计算切线斜率k = f'(x_0)。
3. 根据切线与x 轴相交的方程,求解f(x) = 0,得到新的近似解x_1。
4. 判断x_1 是否满足精度要求,若满足则停止迭代;若不满足,则令x_0 = x_1,返回步骤2。
需要注意的是,牛顿迭代法并不一定能够收敛到方程的根,可能会陷入局部最优解或者发散。
因此,在使用牛顿迭代法时,需要对初始近似解的选择和迭代过程的控制进行合理的调整。
关于数值积分(numerical integration),也称为数值求积,是通过数值计算来求解定积分的方法。
定积分表示曲线与坐标轴之间的面积,常用于求解函数在某个区间上的总体积、质量、电荷等物理量。
数值积分有多种方法,常见的包括梯形法则、辛普森法则、龙贝格法等。
这些方法的基本思想都是将定积分转化为对函数在一系列离散点上的取值进行计算。
以梯形法则为例,其基本思想是将积分区间等分成多个小区间,然后用每个小区间上的函数值构成的梯形的面积来近似表示积分的结果。
具体步骤如下:1. 将积分区间[a, b] 等分成n 个小区间,每个小区间的宽度为h = (b - a) / n。
2. 在每个小区间上计算函数的取值,得到函数在离散点上的值f(x_0), f(x_1), ..., f(x_n)。
数值积分与微分方程数值解法数值积分和微分方程数值解法是数值计算中的重要组成部分,在科学计算、工程分析和实际问题求解中起着不可或缺的作用。
本文将介绍数值积分的基本概念和常用方法,以及微分方程数值解法的应用和实现过程。
一、数值积分的基本概念和常用方法数值积分是求解定积分近似值的方法,通过将连续函数的积分转化为离散形式的求和,以达到近似计算的目的。
常用的数值积分方法包括矩形法、梯形法、辛普森法等。
(1)矩形法:将积分区间等分为若干子区间,然后在每个子区间内取点,用函数在相应点处的取值近似代替该子区间内的函数值,最后将所有子区间的函数值相加得到近似积分值。
(2)梯形法:与矩形法类似,但是将每个子区间近似为一个梯形,通过计算梯形的面积来近似计算积分值。
(3)辛普森法:将积分区间等分为若干子区间,然后在每个子区间内取三个点,根据这三个点构造出一个二次函数,并用该二次函数的积分来近似计算积分值。
二、微分方程数值解法的应用和实现过程微分方程数值解法是对微分方程进行近似求解的方法,通过离散化微分方程来构造数值格式,然后通过数值计算来求解。
常用的微分方程数值解法包括常微分方程的欧拉法、改进欧拉法和龙格-库塔法,以及偏微分方程的有限差分法、有限元法等。
(1)常微分方程数值解法:- 欧拉法:根据微分方程的定义,将微分项近似为差分项,通过迭代逼近真实解。
- 改进欧拉法:在欧拉法的基础上,通过利用两个点的斜率来逼近解的变化率,提高精度。
- 龙格-库塔法:通过多次迭代,根据不同的权重系数计算不同阶数的近似解,提高精度。
(2)偏微分方程数值解法:- 有限差分法:将偏微分方程中的一阶和二阶导数近似为差分项,通过离散化区域和时间来构造矩阵方程组,然后通过求解线性方程组来获得数值解。
- 有限元法:将区域进行剖分,将偏微分方程转化为变分问题,通过选取适当的试函数和加权残差法来逼近真实解。
总结:数值积分和微分方程数值解法是数值计算中重要的工具,能够帮助我们处理实际问题和解决科学工程中的复杂计算。
几种定积分的数值计算方法摘要:本文归纳了定积分近似计算中的几种常用方法,并着重分析了各种数值方法的计算思想,结合实例,对其优劣性作了简要说明.关键词:数值方法;矩形法;梯形法;抛物线法;类矩形;类梯形Several Numerical Methods for Solving Definite Integrals Abstract:Several common methods for solving definite integrals are summarized in this paper. Meantime, the idea for each method is emphatically analyzed. Afterwards, a numerical example is illustrated to show that the advantages and disadvantages of these methods.Keywords:Numerical methods, Rectangle method, Trapezoidal method, Parabolic method, Class rectangle, Class trapezoid1. 引言在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数)(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用. 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数)(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用.另外,对于求导数也有一系列的求导公式和求导法则.但是,在实际问题中遇到求积分的计算,经常会有这样的情况:(1)函数)(x f 的原函数无法用初等函数给出.例如积分 dx e x ⎰-102, ⎰10sin dx xx等,从而无法用牛顿-莱布尼茨公式计算出积分。
积分方程的数值解法及其应用积分方程是一种重要的数学工具,广泛应用于科学和工程等各个领域。
然而,积分方程通常没有解析解,需要借助数值方法来求解。
本文将介绍积分方程的数值解法及其应用。
积分方程的数值解法积分方程的数值解法有很多种,常用的方法包括:•格点法:将积分方程离散化为一组代数方程组,然后用数值方法求解代数方程组。
格点法是积分方程数值解法中最简单的方法,但精度不高。
•边界元法:将积分方程转化为一组边界积分方程,然后用数值方法求解边界积分方程。
边界元法比格点法精度更高,但计算量更大。
•谱法:将积分方程转化为一组谱方程,然后用数值方法求解谱方程。
谱法是一种高精度的积分方程数值解法,但计算量非常大。
积分方程的应用积分方程在科学和工程等各个领域都有广泛的应用,例如:•电磁学:积分方程可以用来求解电磁场问题,如天线设计、微波电路设计等。
•流体力学:积分方程可以用来求解流体力学问题,如流体流动、湍流、热传导等。
•固体力学:积分方程可以用来求解固体力学问题,如弹性力学、塑性力学、断裂力学等。
•化学工程:积分方程可以用来求解化学工程问题,如反应器设计、传质、传热等。
•生物学:积分方程可以用来求解生物学问题,如种群动态、流行病学、药物动力学等。
积分方程数值解法的发展前景积分方程数值解法是一个不断发展的领域,随着计算技术的进步,积分方程数值解法的方法和精度也在不断提高。
近年来,积分方程数值解法在以下几个方面取得了重大进展:•快速算法的开发:近年来,人们开发了许多快速算法来求解积分方程,如快速多极子算法、快速边界元算法、快速谱法等。
这些算法大大提高了积分方程数值解法的速度和效率。
•并行算法的开发:随着并行计算技术的兴起,人们也开发了许多并行算法来求解积分方程。
这些算法可以充分利用多核处理器和分布式计算资源,进一步提高积分方程数值解法的速度和效率。
•自适应算法的开发:自适应算法是一种根据积分方程的局部误差来调整计算精度的算法。
数值计算中的数值积分方法数值计算是应用数学的一个分支,它主要涉及数值计算方法、算法和数值实验。
其中,数值积分作为数值计算中的一个重要环节,其作用在于将连续函数转化为离散的数据,从而方便计算机进行计算和处理。
本文将介绍数值积分的概念、方法和应用。
一、数值积分的概念数值积分是利用数值方法对定积分进行估计的过程。
在数值积分中,积分被近似为离散区间的和,从而可以被计算机进行处理。
数值积分中,被积函数的精确的积分值是无法计算的,而只能通过数值方法进行估计。
数值积分的目的是通过选取合适的算法和参数来尽可能减小误差,达到精度和效率的平衡。
二、数值积分的方法1. 矩形法矩形法是数学上最简单的数值积分方法之一。
矩形法的算法是将要积分的区间分为若干个小区间,然后计算每个小区间中矩形的面积,最后将所有小矩形的面积加起来得到近似的积分值。
矩形法的精度一般较低,适用于计算不需要高精度的函数积分。
2. 梯形法梯形法是数值积分中常用的一种方法,其原理是将区间分为若干个梯形,并计算每个梯形的面积,最后将所有梯形的面积加起来得到近似的积分值。
梯形法的计算精度较高,但其计算量较大。
3. 辛普森法辛普森法是数值积分中一种高精度的方法,它是利用二次多项式去估计原函数。
辛普森法的原理是将区间分为若干等分小区间,并计算每个小区间中的二次多项式的积分值,最后将所有小区间的积分值加起来得到近似的积分值。
辛普森法的优点是其精度高,计算量相对较小。
三、数值积分的应用数值积分方法在各个领域都有广泛的应用。
例如,它可以被用于工程学、物理学和金融学中的数值计算。
在工程学中,数值积分被用于数值模拟和计算机辅助设计中。
在物理学中,数值积分则被用于数值求解微分方程和计算机模拟等领域。
在金融学中,数值积分则被应用于计算复杂的金融模型和风险分析。
总之,数值积分方法是数学和计算机科学中非常重要的一部分。
通过不同的数值积分方法来近似计算定积分,我们能够利用计算机更加高效地进行数学计算和数据分析,从而使得数学和物理等学科的研究者能够更加快速地得出准确的结果。
实验09 数值微积分与方程数值求解(第6章 MATLAB 数值计算)一、实验目的二、实验内容1. 求函数在指定点的数值导数232()123,1,2,3026x x x f x x xx x==2. 用数值方法求定积分(1) 210I π=⎰的近似值。
程序及运行结果:《数学软件》课内实验王平(2) 2221I dx x π=+⎰程序及运行结果:3. 分别用3种不同的数值方法解线性方程组6525494133422139211x y z u x y z u x y z u x y u +-+=-⎧⎪-+-=⎪⎨++-=⎪⎪-+=⎩ 程序及运行结果:4. 求非齐次线性方程组的通解1234123412342736352249472x x x x x x x x x x x x +++=⎧⎪+++=⎨⎪+++=⎩5. 求代数方程的数值解(1) 3x +sin x -e x =0在x 0=1.5附近的根。
程序及运行结果(提示:要用教材中的函数程序line_solution ):(2) 在给定的初值x 0=1,y 0=1,z 0=1下,求方程组的数值解。
23sin ln 70321050y x y z x z x y z ⎧++-=⎪+-+=⎨⎪++-=⎩6. 求函数在指定区间的极值(1) 3cos log ()xx x x xf x e ++=在(0,1)内的最小值。
(2) 33212112122(,)2410f x x x x x x x x =+-+在[0,0]附近的最小值点和最小值。
7. 求微分方程的数值解,并绘制解的曲线2250(0)0'(0)0xd y dyy dx dx y y ⎧-+=⎪⎪⎪=⎨⎪=⎪⎪⎩程序及运行结果(注意:参数中不能取0,用足够小的正数代替):令y 2=y,y 1=y ',将二阶方程转化为一阶方程组:'112'211251(0)0,(0)0y y y x x y y y y ⎧=-⎪⎪=⎨⎪==⎪⎩8. 求微分方程组的数值解,并绘制解的曲线123213312123'''0.51(0)0,(0)1,(0)1y y y y y y y y y y y y =⎧⎪=-⎪⎨=-⎪⎪===⎩程序及运行结果:三、实验提示四、教程:第6章 MATLAB 数值计算(2/2)6.2 数值微积分 p155 6.2.1 数值微分1. 数值差分与差商对任意函数f(x),假设h>0。
数值计算方法数值积分与微分方程数值解数值计算是计算数值结果的一种方法,广泛应用于科学、工程和金融等领域。
数值计算方法涉及到估算数学问题的解,其中包括数值积分和微分方程数值解。
本文将分别介绍数值积分和微分方程数值解的基本原理和常用方法。
一、数值积分数值积分是通过数值计算方法来估计函数的积分值。
积分是数学中的重要概念,广泛应用于物理、经济等领域的问题求解中。
传统的积分计算方法,如牛顿-柯特斯公式和高斯求积法,需要解析求解被积函数,但是对于大多数函数来说,解析求解并不容易或者不可能。
数值计算方法通过离散化被积函数,将积分问题转化为求和问题,从而得到近似的积分结果。
常见的数值积分方法包括梯形法则、辛普森法则和复化求积法。
1. 梯形法则梯形法则是最简单的数值积分方法之一。
它将积分区间划分为若干个小区间,然后在每个小区间上用梯形的面积来近似原函数的面积,最后将所有小区间的梯形面积相加得到近似积分值。
2. 辛普森法则辛普森法则是一种比梯形法则更精确的数值积分方法。
它将积分区间划分为若干个小区间,然后在每个小区间上用一个二次多项式来近似原函数,最后将所有小区间的二次多项式积分值相加得到近似积分值。
3. 复化求积法复化求积法是一种将积分区间进一步细分的数值积分方法。
通过将积分区间划分为更多的小区间,并在每个小区间上应用辛普森法则或者其他数值积分方法,可以得到更精确的积分结果。
二、微分方程数值解微分方程是描述自然现象中变化的数学模型。
求解微分方程的解析方法并不适用于所有的情况,因此需要利用数值计算方法来估计微分方程的解。
常见的微分方程数值解方法包括欧拉法、改进的欧拉法、龙格-库塔法等。
1. 欧拉法欧拉法是最简单的微分方程数值解方法之一。
它通过将微分方程离散化,将微分运算近似为差分运算,从而得到微分方程的近似解。
2. 改进的欧拉法改进的欧拉法是对欧拉法的改进。
它通过使用两个不同的点来估计微分方程的解,从而得到更精确的近似解。
数值计算方法复习知识点数值计算方法是研究计算数值解的方法和数值计算的理论。
它是计算数学的一个分支,主要用于解决无法用解析方法求解的数学模型问题。
本文将综述数值计算方法的一些重要知识点,包括插值与逼近、数值微分与数值积分、线性方程组的直接解法与迭代解法以及常微分方程的数值解法。
一、插值与逼近1.插值:插值是利用已知数据点构造一个函数,使得该函数在给定的数据点上与已知函数完全相等。
常见的插值方法有拉格朗日插值和牛顿插值。
2. 逼近:逼近是从已知数据点构造一个函数,使得该函数在给定的数据点附近与已知函数近似相等。
逼近常用的方法有最小二乘逼近和Chebyshev逼近。
二、数值微分与数值积分1.数值微分:数值微分是通过计算差分商来近似计算函数的导数。
常见的数值微分方法有前向差分、后向差分和中心差分。
2.数值积分:数值积分是通过近似计算定积分的值。
常见的数值积分方法有中矩形法、梯形法和辛普森法。
三、线性方程组的直接解法与迭代解法1.直接解法:直接解法是通过一系列数学运算直接计算线性方程组的解。
常见的直接解法有高斯消元法和LU分解法。
2. 迭代解法:迭代解法是通过迭代计算逼近线性方程组的解的方法。
常见的迭代解法有Jacobi迭代法和Gauss-Seidel迭代法。
四、常微分方程的数值解法1.常微分方程:常微分方程是描述动力系统的数学模型,常用来描述物理系统、生物系统等。
常微分方程的数值解法主要包括初始值问题的一阶常微分方程和常微分方程组的数值解法。
2.常微分方程的数值解法:常微分方程的数值解法有欧拉方法、改进的欧拉方法、龙格-库塔方法等。
这些方法都是将微分方程转化为递推方程,通过迭代计算逼近微分方程的解。
总结:数值计算方法是求解数学模型的重要工具,在科学计算、工程设计和经济管理等领域有广泛的应用。
本文回顾了数值计算方法的一些重要知识点,包括插值与逼近、数值微分与数值积分、线性方程组的直接解法与迭代解法以及常微分方程的数值解法。
(1):方程:取x的范围为[0,10),a=0,b=5;求解步骤:1、对s进行离散,取,将[0,5]中每隔0.01取的数以此记为(i=0,1,2,3….500),即有2、对x进行离散,也取,将[0,10)中每隔0.01取的数以此记为(i=0,1,2,3….999)对于所取的任意有写成矩阵的形式为对于全部的,可写成:1000说明:s的范围可以为x的取值范围,也可以比x的取值范围小,当s的取值范围比x的取值范围小时,可以像上式那样将等式右边第二项的系数矩阵中s取不到的值令k=0,s和x的范围相等仅仅是一种特殊情况。
3、上式可化成F-KF=G的形式,进一步化成AF=G的形式,对其用doolittle分解法求解,可等到一系列离散值。
求解过程中取a=0,b=5,x的范围为[0,10),,,的理论值应为; C++程序如下:#include <iostream.h>#include <fstream.h>#include <math.h>doublea[1000][1000],b[1000][1000],u[1000][1000],l[ 1000][1000],y[1000],f[1000];double m=2.0;double ds=0.01;double si,xj,k,g,xo,gg;double kx(double s,double x){k=s*x*x;return k;}double gx(double x){g=x*x*(-103.1667)-4*x+5;return g;}void main(){fstream outfile;outfile.open("jifenshi.dat",ios::out);for (int j=0;j<1000;j++){xj=j;for (int i=0;i<=500;i++){si=i;a[j][i]=m*ds*kx(si/100,xj/100);}}for (int jj=0;jj<1000;jj++){for (int ii=0;ii<1000;ii++){if (ii==jj)b[jj][ii]=1-a[jj][ii];elseb[jj][ii]=-a[jj][ii];}}for (int c=0;c<1000;c++) //doolittle {u[0][c]=b[0][c];l[c][0]=b[c][0]/u[0][0];}for (int d=1;d<1000;d++){for (int e=d;e<1000;e++){double ff=0;for (int f=0;f<=d-1;f++){ff+=l[d][f]*u[f][e];}u[d][e]=b[d][e]-ff;}for (int h=d;h<1000;h++){double nn=0;for (int n=0;n<=d-1;n++){nn+=l[h][n]*u[n][d];}l[h][d]=(b[h][d]-nn)/u[d][d];}}y[0]=gx(0);for (int o=1;o<1000;o++){xo=o;gg=gx(xo/100);double pp=0;for (int p=0;p<=o-1;p++){pp+=l[o][p]*y[p];}y[o]=gg-pp;}f[999]=y[999]/u[999][999];for (int q=998;q>=0;q--){double tt=0;for (int t=q+1;t<1000;t++){tt+=u[q][t]*f[t];}f[q]=(y[q]-tt)/u[q][q];}for (int z=0;z<1000;z++){outfile<<f[z]<<endl;}outfile.close();}得出的的一系列离散值与理论值在matlab中成图如下:由图可看出两者相差比较小数值计算所得的值与理论值的相对误差图分析:离散的时候步长取的越小,相对误差越小,另外,我在一开始取定和模型时由此算出函数,此函数系数为小数,在取的时候进行了一定的近似处理,因此在由和算时产生了一定的人为误差。
实验七积分方程的数值计算方法
1. 实验描述
计算
32
sin(4)x
x e dx
-
⎰定积分的近似值,起始容差00.00001
ε=
1.用组合梯形公式M=10计算。
2.用组合辛普生公式M=5计算。
3.用龙贝格积分计算。
4.用自适应积分方法计算。
2. 实验内容
1. 用组合梯形公式M=10计算。
图1. 组合梯形算法流程图
将积分区间[],a b划分为宽度为()/
=-的M个子区间,再将各区间的面积
h b a M
求和即可得到近似面积。
2.用组合辛普生公式M=5计算。
Array
图2. 组合辛普森算法流程图
将积分区间[],a b划分为宽度为()/2
h b a M
=-的2M个子区间,再由辛普森公式将各区间的面积求和即可得到近似面积。
3.用龙贝格积分计算。
图3.龙贝格积分算法流程图
①.由递归梯形公式序列得到递归辛普森序列序列。
②.由递归辛普森序列序列得到递归布尔公式序列。
③.通过理查森改进提高误差项的阶数。
4.用自适应积分方法计算
图4.自适应积分算法流程图
在辛普森公式基础上,将区间再进行划分,即为自适应积分。
3. 实验结果及分析
真实值S = 0.19971466216144
1. 用组合梯形公式M=10计算。
面积近似值S1 =0.16965032127666。
误差error1=0.03006434088479。
2.用组合辛普森公式M=5计算
4.660686426147481*10-。
面积近似值S2 =0.19966805529718。
误差error2=5 3.用龙贝格积分计算。
表1.龙贝格积分表
面积近似值S3 =0.19971378242654。
误差error3=5
2.254986251026825*10-。
4.用自适应积分方法计算
面积近似值S4 =0.19971391278871。
误差error4=6
3.238553837777504*10-。
4. 结论
自适应计分方法和龙贝格积分法都可以得到较高的精确度。
组合梯形和组合辛普森可以增大M值,得到较高的精确度。
附件(代码)
1.
f=inline('sin(4*x)*exp(-2*x)','x');%定义函数
A=0;t=0;S=0;
for k=1:9
t=k*0.3;
A=A+0.3*feval(f,t);
end
S=A+0.15*(feval(f,0)+feval(f,3))%组合梯形的近似值
RS=int(sin(4*x)*exp(-2*x),0,3)%真实值
error=RS-S%真实值与近似值的误差
2.
f=inline('sin(4*x)*exp(-2*x)','x');%定义函数
S=0;X=zeros(1,11);h=0.3;
for j=1:11
X(j)=(j-1)*h;
end
for k=1:5
S=S+h/3*(feval(f,X(2*k-1))+4*feval(f,X(2*k))+feval(f,X(2*k+1))); end
RS=int(sin(4*x)*exp(-2*x),0,3);%真实值
S
error=RS-S%真实值与近似值的误差
3.
function [R,quad,err,h]=romber(f,a,b,n,tol)
%INPUT-f is the intergrand input as astring 'f'
% -a and b are upper and lower limits of interation
% -n is the maxium number of rows in the table
% -tol is the tolerance
%OUTPT-R is the Romberg table
% -quad is the quadrature value
% -err is the error estimate
% -h is the smallest step size uesd
M=1;
h=b-a;
err=1;
J=0;
R=zeros(4,4);
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
while((err>tol)&(J<n))|(J<4)
J=J+1;
h=h/2;
s=0;
for p=1:M
x=a+h*(2*p-1);
s=s+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*s;
M=2*M;
for K=1:J
R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);
end
err=abs(R(J,J)-R(J+1,K+1));
end
quad=R(J+1,J+1);
4.
function Z=srule(f,a0,b0,tol0)
%INPUT-f is the intergrand input as astring 'f'
% -a0 and b0 are upper and lower limits of intergration % -tol0 is the tolerance
%OUTPUT-Z is a 1x6vector [a0 b0 S S2 err tol1]
h=(b0-a0)/2;
C=zeros(1,3);
C=feval(f,[a0 (a0+b0)/2 b0]);
S=h*(C(1)+4*C(2)+C(3))/3;
S2=S;
tol1=tol0;
err=tol0;
Z=[a0 b0 S S2 err tol1];
function [SRmat,quad,err]=adapt(f,a,b,tol)
%INPUT-f is the intergrand input as astring 'f'
% -a and b are upper and lower limits of interation % -tol is the tolerance
%OUTPT-SRmat is the Romberg table
% -quad is the quadrature value
% -err is the error estimate
% -Initialize values
SRmat=zeros(30,6);
iterating=0;
done=1;
SRvec=zeros(1,6);
SRvec=srule(f,a,b,tol);
SRmat(1,1:6)=SRvec;
m=1;
state=iterating;
while(state==iterating)
n=m;
for j=n:-1:1
p=j;
SR0vec=SRmat(p,:);
err=SR0vec(5);
tol=SR0vec(6);
if (tol<=err)
%Bisect interval,apply Simpson'rule
%recursively,and determine error
state=done;
SR1vec=SR0vec;
SR2vec=SR0vec;
a=SR0vec(1);
b=SR0vec(2);
c=(a+b)/2;
err=SR0vec(5);
tol=SR0vec(6);
tol2=tol/2;
SR1vec=srule(f,a,c,tol2);
SR2vec=srule(f,c,b,tol2);
err=abs(SR0vec(3)-SR1vec(3)-SR2vec(3))/10;
%Accuracy test
if (err<tol)
SRmat(p,:)=SR0vec;
SRmat(p,4)=SR1vec(3)+SR2vec(3);
SRmat(p,5)=err;
else
SRmat(p+1:m+1,:)=SRmat(p:m,:);
m=m+1;
SRmat(p,:)=SR1vec;
SRmat(p+1,:)=SR2vec;
state=iterating;
end
end
end
end
quad=sum(SRmat(:,4));
err=sum(abs(SRmat(:,5)));
SRmat=SRmat(1:m,1:6);。