变步长复化辛普森公式计算积分
0 2. 编写用变步长复化辛普森公式计算积分
b
f (x)dx 的程序。
a
用上面编写的程序计算下列积分并分析计算结果
（1 ）
x cos xdx （ 2 ） 0 2 x 2 cosx 2dx （ 3） 1 xdx
程序：
function S=bianfuhuasimpson(fx,a,b,eps,M)
% 变步长复合 simpson 求积公式
% 调用方式： S=fuhuasimpson(@fx,a,b,epsilon)
% fx -- 求积函数（函数文件）
% a, b -- 求积区间
% eps -- 计算精度
% M-- 最大允许输出划分数
n=1;
h=(b-a)/n;
T1=h*(feval(fx,a)-feval(fx,b))/2;
Hn=h*feval(fx,(a+b)/2);
S1=(T1+2*Hn)/3;
n=2*n;
% 最好与倒数第三行保持一致（变步长）
while n<=M
T2=(T1+Hn)/2;
Hn=0;
h=(b-a)/n;
for j=1:n
x(j)=a+(j-1/2)*h;
y(j)=feval(fx,x(j));
Hn=Hn+y(j);
end
Hn=h*Hn;
S2=(T2+2*Hn)/3;
fprintf(' n=%2d S2=%-12.9f S2-S1=%-12.9f\n',n,S2,abs(S2-S1)); if abs(S2-S1) break; else T1=T2; S1=S2; n=2*n; end end S=S2; % 达到下列条件之一，则运算终止: % (1).abs(S2-S1) % (2). 下一次的n>M % 输入1 ：S=bianfuhuasimpson(inline('sqrt(x)*cos(x)'),0,pi,10e-6,2000) % 输入2 ：S=bianfuhuasimpson(inline('2*x^2*cos(x^2)'),0,sqrt(pi),10e-6,2000) % 输入3 ：S=bianfuhuasimpson(inline('sqrt(x)'),0,1,10e-6,2000) 输出结果： （1 ） S=bianfuhuasimpson(inline('sqrt(x)*cos(x)'),0,pi,10e-6,2000) n= 2 S2=-0.016369112 S2-S1=0.944423778 n= 4 S2=-0.450266122 S2-S1=0.433897010 n= 8 S2=-0.669839370 S2-S1=0.219573248 n=16 S2=-0.781318443 S2-S1=0.111479074 n=32 S2=-0.837710689 S2-S1=0.056392245 n=64 S2=-0.866141900 S2-S1=0.028431211 n=128 S2=-0.880440980 S2-S1=0.014299080 n=256 S2=-0.887620063 S2-S1=0.007179083 n=512 S2=-0.891220052 S2-S1=0.003599989 n=1024 S2=-0.893023740 S2-S1=0.001803689 S = -0.8930 （2 ） S=bianfuhuasimpson(inline('2*x^2*cos(x^2)'),0,sqrt(pi),10e-6,2000) n= 2 S2=1.076354541 S2-S1=2.092222287 n= 4 S2=0.039359358 S2-S1=1.036995183 n= 8 S2=-0.430456535 S2-S1=0.469815894 n=16 S2=-0.662796649 S2-S1=0.232340113 n=32 S2=-0.778823323 S2-S1=0.116026674 n=64 S2=-0.836827971 S2-S1=0.058004648 n=128 S2=-0.865829756 S2-S1=0.029001785 n=256 S2=-0.880330615 S2-S1=0.014500859 n=512 S2=-0.887581042 S2-S1=0.007250427 n=1024 S2=-0.891206256 S2-S1=0.003625214 S = -0.8912 （3 ） S=bianfuhuasimpson(inline('sqrt(x)'),0,1,10e-6,2000) n= 2 S2=0.489859598 S2-S1=0.185121744 n= 4 S2=0.579745947 S2-S1=0.089886349 n= 8 S2=0.623731522 S2-S1=0.043985575 n=16 S2=0.645384849 S2-S1=0.021653327 n=32 S2=0.656091436 S2-S1=0.010706587 n=64 S2=0.661402273 S2-S1=0.005310836 n=128 S2=0.664042680 S2-S1=0.002640407 n=256 S2=0.665357576 S2-S1=0.001314896 n=512 S2=0.666013147 S2-S1=0.000655572 n=1024 S2=0.666340270 S2-S1=0.000327122 S = 0.6663 Welcome To Download !!! 欢迎您的下载，资料仅供参考！ 实验四、复化梯形公式和复化Simpson公式的精度比较 （2学时） 一、实验目的与要求 1、熟悉复化Simpson公式和复化梯形公式的构造原理； 2、熟悉并掌握二者的余项表达式； 3、分别求出准确值，复化梯形的近似值，复化Simpson的近似值，并比较后两 者的精度； 4、从余项表达式，即误差曲线，来观察二者的精度，看哪个更接近于准确值。 二、实验内容： 对于函数 sin () x f x x =，试利用下表计算积分1 sin x I dx x =?。 表格如下： 注：分别利用复化梯形公式和复化Simpson公式计算，比较哪个精度更好。其中：积分的准确值0.9460831 I=。 三、实验步骤 1、熟悉理论知识，并编写相应的程序； 2、上机操作，从误差图形上观察误差，并与准确值相比较，看哪个精度更好； 3、得出结论，并整理实验报告。 四、实验注意事项 1、复化梯形公式，程序主体部分： for n=2:10 T(n)=0.5*T(n-1) for i=1:2^(n-2) T(n)=T(n)+(sin((2*i-1)/2^(n-1))/((2*i-1)/2^(n-1)))/2^(n-1)； end end 2、复化Simpson公式，程序主体部分： for i=1:10 n=2.^i x=0:1/n:1 f=sin(x)./x f(1)=1 s=0 for j=1:n/2 s=s+f(2*j) end t=0 for j=1:(n/2-1) t=t+f(2*j-1) end S(i)=1/3/n*(f(1)+4*s+2*t+f(n+1)) end 五．实验内容 复化梯形公式和复化辛普森公式的引入 复化梯形公式： 1 10[(()]2 n n k k k h T f x f x -+==+∑； 复化辛普森公式： 1 1102 [(4()()]6n n k k k k h S f x f x f x -++ ==++∑； 根据题意和复化梯形公式、复化辛普森公式的原理编辑程序求解代码如下： Matlab 代码 clc s=quad('sin(x)./x',0,1) p1=zeros(10,1); 实验 目 的 或 要 求1、利用复化梯形公式、复化simpson 公式计算积分 2、比较计算误差与实际误差 实 验 原 理 （ 算 法 流 程 图 或 者 含 注 释 的 源 代 码 ） 取n=2,3,…,10分别利用复化梯形公式、复化simpson 公式计算积分1 20I x dx =?，并与真值进行比较，并画出计算误差与实际误差之间的曲线。 利用复化梯形公式的程序代码如下： function f=fx(x) f=x.^2; ％首先建立被积函数，以便于计算真实值。 a=0; ％积分下线 b=1; ％积分上线 T=[]; ％用来装不同n 值所计算出的结果 for n=2:10; h=(b-a)/n; ％步长 x=zeros(1,n+1); ％给节点定初值 for i=1:n+1 x(i)=a+(i-1)*h; ％给节点赋值 end y=x.^2; ％给相应节点处的函数值赋值 t=0; for i=1:n t=t+h/2*(y(i)+y(i+1)); ％利用复化梯形公式求值 end T=[T,t]; ％把不同n 值所计算出的结果装入 T 中 end R=ones(1,9)*(-(b-a)/12*h.^ 2*2); ％积分余项（计算误差） true=quad(@fx,0,1); ％积分的真实值 A=T-true; ％计算的值与真实值之差（实际误差） x=linspace(0,1,9); plot(x,A,'r',x,R,'*') ％将计算误差与实际误差用图像画出来 注：由于被积函数是x.^2，它的二阶倒数为2，所以积分余项为：(-(b-a)/12*h.^ 2*2) 实 验 原 理 （ 算 法 流 程 图 或 者 含 注 释 的 源 代 码）利用复化simpson 公式的程序代码如下： 同样首先建立被积函数的函数文件： function f=fx1(x) f=x.^4; a=0; ％积分下线 b=1; ％积分上线 T=[]; ％用来装不同n值所计算出的结果 for n=2:10 h=(b-a)/(2*n); ％步长 x=zeros(1,2*n+1); ％给节点定初值 for i=1:2*n+1 x(i)=a+(i-1)*h; ％给节点赋值 end y=x.^4; ％给相应节点处的函数值赋值 t=0; for i=1:n t=t+h/3*(y(2*i-1)+4*y(2*i)+y(2*i+1)); ％利用复化simpson公式求值end T=[T,t] ; ％把不同n值所计算出的结果装入T中 end R=ones(1,9)*(-(b-a)/180*((b-a)/2).^4*24) ; ％积分余项（计算误差） true=quad(@fx1,0,1); ％积分的真实值 A=T-true; ％计算的值与真实值之差（实际误差） x=linspace(0,1,9); plot(x,A,'r',x,R,'*') （一） 实验目的 熟悉并掌握数值积分的方法，重要训练复化梯形公式，复化simpson 公式以及romberg 积分。 （二） 问题描述 问题三数值积分椭圆周长的计算。考虑椭圆22221x y a b +=,为计算其周长,只要计算其第一象限的长度即可. 用参数方程可以表示为cos (0/2)sin x a t t y b t π=?≤≤?=? , 计算公式为/0π? 为计算方便,我们可以令1a =,即计算下面的积分 / 0π?/0π=? (/0π?/0a π=?可以归结为上面的形式) 采用复化梯形公式，复化Simpson 公式以及Romberg 积分的方法计算积分 / 0()I b π=? 给出通用程序，该通用程序可以计算任何一个函数在任意一个区间在给定的精度下的数值积分。程序输出为计算出的数值积分值以及计算函数值的次数。 （三） 算法介绍 首先利用给出的各迭代公式，设计程序。在matlab 对话框中输入要计算的函数，给出区间和精度。 复化梯形的迭代公式为： ; 复化simpson迭代公式为： ; Romberg迭代公式为： 。 （四）程序 对于复化梯形公式和复化simpson公式，我们放在中。 （%标记后的程序可用来把b看为变量时的算法实现） %复化梯形公式 function y=jifenn(f,n,a,b) （说明：f表示任一函数，n精度，a，b为区间）fi=f(a)+f(b); h=(b-a)/n; d=1; %function f=jifen(n,a,b,c) %syms t %y=sqrt(1+(c^2-1)*cos(t)^2); %ya=subs(y,t,a); %yb=subs(y,t,b); %fi=ya+yb; for i=1:n-1 x=a+i*h; fi=fi+2*f(x); 龙格库塔法 龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。 龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。 龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。其中4 阶龙格库塔法是最常用的一种方法。因为它相当精确、稳定、容易编程。在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。如果需要较高的精度, 可采取减小步长的方法即可。4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。 1、初值问题 对于一阶常微分方程的初值问题 根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。 2、离散化 取步长h=(b-a)/n，将区间[a , b]分成n个子区间： a=<=b 在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲 线上至少有一点，它的斜率与整段曲线的平均斜率相同， 得=y’() (0<<1) 其中，= 可以将上式改写成y()=y()+h*K (2.1) 其中K为平均斜率，K=f() 公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。 欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。 3、Euler法 欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙 格库塔法的基础。 首先，令、为y() 及y()的近似值，并且令平均斜 率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1） 4、改进的欧拉法 此种方法是取、两点的斜率的平均值作为平均斜率K， 即K= ,其中、均为y()以及y()的近似值，就得到 改进后的欧拉公式（4.1） 数值计算方法上机题目3 一、计算定积分的近似值： 2 21x e xe dx =? 要求： （1）若用复化梯形公式和复化Simpson 公式计算，要求误差限7102 1 -?=ε，分别利用他们的余项估计对每种算法做出步长的事前估计； （2）分别利用复化梯形公式和复化Simpson 公式计算定积分； （3）将计算结果与精确解比较，并比较两种算法的计算量。 1.复化梯形公式 程序： 程序1（求f （x ）的n 阶导数： syms x f=x*exp(x) %定义函数f （x ） n=input('输入所求导数阶数:') f2=diff(f,x,n) %求f(x)的n 阶导数 结果1 输入n=2 f2 = 2*exp(x) + x*exp(x) 程序2： clc clear syms x%定义自变量x f=inline('x*exp(x)','x') %定义函数f(x)=x*exp(x)，换函数时只需换该函数表达式即可 f2=inline('(2*exp(x) + x*exp(x))','x') %定义f(x)的二阶导数,输入程序1里求出的f2即可。 f3='-(2*exp(x) + x*exp(x))'%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值 e=5*10^(-8) %精度要求值 a=1 %积分下限 b=2 %积分上限 x1=fminbnd(f3,1,2) %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值 for n=2:1000000 %求等分数n Rn=-(b-a)/12*((b-a)/n)^2*f2(x1) %计算余项 if abs(Rn) syms t; h=(b-a)/n; f=0;l=0; for k=1:n+1 x(k)=a+(k-1)*h w(k)=0.5*h+x(k) z(k)=subs(y,t,x(k)); end for i=2:n f=f+z(i); end for i=1:n l=l+w(i); end q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T=h/2*(q+p+2*f); T=vpa(T,7) clc,clear; syms t; a=1;b=2; y=exp(1/t); n=5; simpson(a,b,n,y); （3）变步长求积法：以书本例4.5为例function [f]=TN(a,b,y,R0) syms t; T=[]; f=0; q=subs(y,t,a); if y=='sin(t)/t'&&a==0 q=1; end p=subs(y,t,b); T(1)=(b-a)/2*(q+p); i=2; n=i-1; h=(b-a)/n; z1=a+h/2; z2=subs(y,t,z1); 摘要 在工程实验及研究中，实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系.可以说，曲线拟合模型与我们的生活生产密切相关. 本课题着重介绍曲线拟合模型及其应用，其中包括它的基本思想、模型的建立、以及具体应用.为了更好的了解曲线拟合模型，可以将它分为线性与非线性模型，在模型建立的基础上我们可以用最小二乘法来解决一些我们日常所应用的问题. 关键词曲线拟合；线性与非线性模型；最小二乘发 目录 引言 (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) 实验五 一、实验名称 复合梯形积分和复合Simpson 积分计算数值积分 二、实验目的与要求： 实验目的: 掌握复合梯形积分和复合Simpson 积分算法。 实验要求：1.给出复合梯形积分和复合Simpson 积分算法思路， 2.用C 语言实现算法，运行环境为Microsoft Visual C++。 三、算法思路： 我们把整个积分区间[a,b]分成n 个子区间[xi,xi+1]，i=0，1,2,…,n，其中x0=a ，xn+1=b 。这样求定积分问题就分解为求和问题： ?∑?=-==b a n i x x i i dx x f dx x f S 11)()( 当这n+1个结点为等距结点时，即n a b h ih a x i /)(-=+=，其中，i=0,1,2,…,n ，复化梯形公式的形式是 ∑=-+=n i i i n x f x f h S 1 1)]()([2 算法： input n 0.0←S for i=1 to n do ))()((2 1i i x f x f h S S ++ ←- end do output S 如果n 还是一个偶数，则复合Simpson 积分的形式是 ∑=--++=2 /1 21222)]()(4)([3n i i i i n x f x f x f h S 算法： input n 0.0←S for i=1 to n/2 do ))()(4)((3 21222i i i x f x f x f h S S +++ ←-- end do output S 四、实验题目： 五、问题的解： 编写程序（程序见后面附录），输出结果如下： 关于复化辛普森(simpson)公式在线路坐标计算中的应用 天津西站项目部刘思传 摘要：本文里利用辛普森公式导证了线路坐标计算的公式，并在卡西欧FX-4800P计算器中编写了中边线坐标计算的源程序。 关键词：复化辛普森公式，线路坐标计算，曲率。 一.引言 随着我国道路建设等级和质量水平的飞速发展，公路、铁路建设的机械化和日产量日益提高，促使施工中在满足设计精度的前提下，尽可能快速、准确地进行测量放样和检查工作，本文线路曲率变化的特点，利用复化辛普森公式导证了线路坐标计算的通用公式，并利用卡西欧FX-4800P计算器编写了计算线路中边线坐标的源程序。 二.复化辛普森公式数学模型 把积分区间分成偶数等分，记，其中是节点总数，是积分子区间的总数。 记，，在每个区间上用辛普森数值积分公式计算，则得到复化辛普森公式，记为。 复化辛普森积分计算公式 而，称 (1) 式(1)即为辛普森复化公式。 三.线路坐标计算 2. 回旋曲线上点位坐标方位角的计算 如图1，设回旋曲线起点A 的曲率为A ρ，其里程为DK A ；回旋曲线终点B 的曲率为B ρ，其里程为DK B ，Ax ’'y 为以A 为坐标原点，以A 点切线为'x 轴的局部坐标系；Axy 为线路坐标系。 由此回旋曲线上各点曲率半径为R i 和该点离曲线起点的距离?i 成反比，故此任意点的曲率为 c l R i i i /1==ρ（=为常数）. （2） y ' Y B 图1 由式（2）可知，回旋曲线任意点的曲率按线性变化，由此回旋曲线上里程为DK i 点的曲率为 )(A i A B A B A i DK DK DK DK ---+=ρρρρ （3） 当曲线右偏时，取正；当曲线左偏时取负。在图1中有 ???????=== ?I A DK DK i i i dl dl dl R d ρβρβ1 （4） 将式（3）代入式（4）得 πρρβ180 *)(2A i A i i DK DK -+= （5） 若已知回旋曲线起点A 在线路坐标系下切线坐标方位角αA ，则里程为Dk i 点切线坐标方位角为 i A i βαα+= π180 （6） 将式(5)代入式（6）得 *)(2A i A i A i DK DK -++=ρραα π180 （7） 对于式(7) ，当，时，，则a i =a A ，式(7)变成计算直线段上任意点切线坐标方位角计算公式；当，时，， ，则式(7)代表圆曲线上任意点切线坐标方位角 计算公式。 可见，若已知曲线段起点和终点的曲率及起点的切线坐标方位角，式(7)便能计算任意线型点位切线坐标方位角。 3、回旋曲线点位坐标计算 由图1可得回旋曲线上点位在坐标系下坐标计算公式： 数值分析第五次程序作业 PB09001057 孙琪 【问题】 分别编写用复化Simpson积分公式和复化梯形积分公式计算积分的通用程序；用如上程序计算积分： 取节点并分析误差； 简单分析你得到的数据。 【复化Simpson积分公式】 Simpson法则： 使用偶数个子区间上的复合Simpson法则： 设n是偶数， 则有 将Simpson法则应用于每一个区间，得到复合Simpson法则： 公式的误差项为： 其中δ 【复化梯形积分公式】 梯形法则：对两个节点相应的积分法则称为梯形法则： 如果划分区间[a,b]为： 那么在每个区间上可应用梯形法则，此时节点未必是等距的，由此得到复合梯形法则： 对等间距h=(b-a)/n及节点，复合梯形法则具有形式： 误差项为： 【算法分析】 复合Simpson法则和复合梯形法则的算法上述描述中都已介绍了，在此不多做叙述。 【实验】 通过Mathematica编写程序得到如下结果： 利用复化Simpson积分公式得： 可以看出，当节点数选取越来越多时，误差项越来越小，这从复合的Simpson公式很好看出来，因为在每一段小区间内，都是用Simpson法则去逼近，而每一段的误差都是由函数在该区间内4阶导数值和区间长度的4次方乘积决定的，当每一段小区间越来越小时，相应的每一段小区间内的逼近就会越来越好，从而整体的逼近效果就会越来越好。 利用复化梯形积分公式得： 可以看出，当节点数选取越来越多时，误差项越来越小，这从复合的梯形公式很好看出来，因为在每一段小区间内，都是用梯形法则去逼近，而每一段的误差都是由函数在该区间内2阶导数值和区间长度的2次方乘积决定的，当每一段小区间越来越小时，相应的每一段小区间内的逼近就会越来越好，从而整体的逼近效果就会越来越好。 【分析】 通过对上述两种法则的效果来看，复合Simpson法则的误差要比复合梯形法则收敛到0更快，说明复合Simpson法则逼近到原来的解更快，这主要是因为在每一段小区间内，复合Simpson法则利用得是Simpson法则，复合梯形法则利用得是梯形法则，前者的误差项要比后者的误差项小很多，因此造成了逼近速度的不一样。 第五章 常微分方程的数值解法 5.2 龙格-库塔法 一、教学目标及基本要求 通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。 二、教学内容及学时分配 本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性 三、教学重点难点 1．教学重点：龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。 四、教学中应注意的问题 多媒体课堂教学为主。适当提问，加深学生对概念的理解。 五、正文 龙格－库塔方法 引言 龙格－库塔方法的基本思想 '11()() ()()()(,())n n n n y x y x y y x y x hf y h ξξξ++-=?=+ 令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。 欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。 考察改进的欧拉公式 11211 ()22 n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++ 它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉公式预报得到。 可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间 1[,]n n x x +的平均斜率*K 。这就是龙格－库塔方法的基本思想。 1、二阶龙格－库塔方法 考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率 12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+ 取1(,)n n K f x y =，如何预报n p x +处斜率2K ？ 仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +： 1n p n y y phK +=+ 然后用n p y +计算2(,)n p n p K f x y ++=，从而得 112121[(1)](,)(,) n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+ 适当选取,p λ，使上述公式具有较高得精度。假定()n n y y x =，分别将12 ,K K 泰勒展开： '1(,)()n n n K f x y y x == 212 ' '' 2 (,) (,)[(,)(,)(,)]()()()() n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++ 代入得 '2''31()()()()n n n n y y x hy x ph y x O h λ+=+++ 按泰勒展开2' '' 31()()()()2 n n n n h y y x hy x y x O h +=+++ 比较得，只要1 2 p λ= ，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==， 12n n y y hk +=+，1(,)n n k f x y =，21(,)22 n n h h k f x y k =++ Simpson算法及其推广形式 摘要：本文研究了辛普森公式的数值积分的计算方法问题，并且更进一步研究了变步长复化的辛普森公式和二重积分的辛普森公式的问题。首先是对 一维辛普森公式和变步长复化辛普森公式以及二维辛普森公式的推导及 其算法，进行误差分析，并且列举了实例。然后，对辛普森公式进行改 进，这里的改进最主要是对辛普森公式的代数精度进行提高，从而使辛 普森公式对积分的计算更加精确。另外，还研究了辛普森公式的推广形 式。最后，在结论的当中列举了一个例子。 关键词：辛普森公式算法改进推广形式二重积分的辛普森公式 Abstract：This paper first studies the calculation methods of the numerical integration in simpson formula, and then study of the long-simpson formula and the double integral simpson formula problem. First, study the algorithm and derived of one-dimensional simpson formula and step-change in simpson formula, as well as two-dimensional simpson formula, and then analysis the error. Finally , list the example. In this , improve the simpson formula. This improved the most important is to incre ase the simpson formula’s accuracy of algebra. Besides, we study the simpson formula’s promotion of forms. At the last, we list a example in the conclusion. Key word:The simpson formula, Algorithm, Improve, Promotion of forms, The simpson formula of the two-dimensional integral. 数值分析实验 三 班级：10信计2班 学号：59 姓名：王志桃 分数 一·问题提出： 选用复合梯形公式，复合Simpson 公式，计算 （1） I =dx x ?-4 10 2sin 4 ()5343916.1≈I （2） I = dx x x ?1 sin ()9460831.0,1)0(≈=I f （3） I = dx x e x ?+1 024 （4） I = () dx x x ?++1 021 1ln 二·实验要求： 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 复化梯形公式与辛普森公式的matlab程序【程序代码】 cclc; disp('1.复化梯形公式求解'); disp('2.simpson公式求解'); disp('请进行选择：'); c=input(' '); if c==1 clc; disp('复化梯形公式'); disp('请输入积分下限'); a=input('a='); disp('请输入积分上限'); b=input('b='); disp('请输入等分的数目'); n=input('n='); h=(b-a)/n; s1=0; for i=1:n-1 s1=s1+fun1(i*h); end disp('复化梯形公式的结果：'); T=h/2*(fun1(a)+2*s1+fun1(b)) else if c==2 clc; disp('simpson公式'); disp('请输入积分下限'); a=input('a='); disp('请输入积分上限'); b=input('b='); disp('请输入等分的数目'); n=input('n='); h=(b-a)/n; s2=0; for i=0:n-1 s2=s2+fun1((i+0.5)*h); end disp('辛普森公式的结果：'); S=h/6*(fun1(a)+4*s2+2*s1+fun1(b)) end end disp('菜单选项'); disp('1.继续运算'); disp('2.退出程序！'); p=input(' '); if p==1 (fuhua); else if p==2 disp('正在退出，请稍候。。。'); 《复合梯形公式》实验报告 实验名称：验证数值积分求积公式及复合梯形公式程序设计成绩：___________ 专业班级：数学与应用数学1202班姓名：王晓阳学号：2012254010228 实验日期：2014 年10月20日 实验报告日期：2014年11月3日 一．实验目的 1掌握定积分的数值求解方法，验证数值积分求积公式. 2.掌握数值积分的基本思想. 3.掌握matlab实现数值积分函数的调用格式. 4.编写复合梯形公式matlab程序及学会调用. 5.学会用复合梯形公式求函数近似解. 二、实验内容 1.数值积分的实现 (1)被积函数是一个解析式 Matlab提供了quad函数和quadl函数来求定积分.它们的调用格式为： Quad(filename,a,b,tol,trace) Quadl(filename,a,b,tol,trace) 其中filename 是被积函数名。a 和b 分别是定积分的下限和上限。Tol 用来控制积分精度，默认时取610tol -=。Trace 控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0。 例6.20 用两种不同的方法求2 10x I e dx -=?. (2)被积函数由一个表格定义 在matlab 中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X 、Y 定义函数关系Y=f(X).X 、Y 是两个等长的向量； ()12,,n X x x x =，()12,,n Y y y y =，并且12n x x x <<<，积分区间是[]1,n x x 。 例6.21用trapz 函数计算210 x I e dx -=?. (3)二重积分数值求解 Matlab 提供的dblquad 函数可以直接求出二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace) 该函数求f(x,y)在[a,b]*[c,d]区域上的二重定积分。 例6.22 计算二重定积分()2 122212sin x e x y dxdy ---+? ? 2.复合梯形公式 由于牛顿-柯特斯公式在8n ≥时不具有稳定性，故不可能通过提高阶的方法来提高求积精度.为了提高精度通常可把积分区间分成若干子区间(通常是等分)，再在每个子区间上用低阶求积公式，这种方法称为复合求积法. 将区间[],a b 划分为n 等份，分点0,0,1, k x x kh k n =+=，0,,n b a a x b x h n -=== ，在每个子区间[]1,k k x x +上采用梯形公式 四阶龙格一库塔法 通常所说的龙格一库塔法是指四阶而言的．我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格一库塔公式 (9.22) 公式(9.22)的局部截断误差的阶为． 龙格一库塔法具有精度高，收敛，稳定(在一定的条件下)，计算过程中可以改变步长，不需要计算高阶导数值等优点．但仍需计算在一些点上的值，如四阶龙格－库塔法每计算一步需要算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”．(即开始几点的近似值)．例3.用标准龙格一库塔法求初值问题 在处的解． 解因与．若应用标准龙格一库塔方法公式(9.22)计算，对于n=0时，则有 于是得 这个值与准确解在处的值已十分接近．再对n=1，2，3，4应用式(9.22)计算，具体计算结果如表3所示： 例3写出用四阶龙格――库塔法求解初值问题 的计算公式，取步长h＝0.2计算y(0.4)的近似值。至少保留四位小数。 解此处f(x，y)＝8－3y，四阶龙格――库塔法公式为 其中κ1＝f(x k，y k)；κ2＝f(x k+0.5h，y k+0.5hκ1)； κ3＝f(x k+0.5h，y k+0.5hκ2)；κ4＝f(x k+h，y k+hκ3) 本例计算公式为： 其中κ1＝8－3y k；κ2＝5.6－2.1y k； κ3＝6.32－2.37y k；κ4＝4.208－1.578y k ＝1.2016＋0.5494y k (k＝0，1，2，…) 当x0＝0，y0＝2， y(0.2)≈y1＝1.2016＋0.5494y0＝1.2016＋0.5494×2＝2.3004 y(0.4)≈y2＝1.2016＋0.5494y1＝1.2016＋0.5494×2.3004＝2.4654 2012-2013（1）专业课程实践论文 复化梯形公式 余良，0818180125，R数学08-1班 一、算法理论 求?baf ()?x baf 的值是科学技术中经常遇到的计算问题 ,但在很多情况下 ,()x f 的原函数不易求得 ,或非常复杂 ;另外 ,在一些应用中 ,函数 ()x f 是用函数表形式给出而没有解析式 ,这时就采用定积分的数值计算方法 ,以解决定积分的近似计算。数值求积公式的应用不仅在近似计算本身 ,在初等数学中 ,某些数列求前 n 项之和公式的推导颇为繁琐 ,应用复化梯形公式可方便地导出这些公式。复化梯形公式用牛顿——莱不尼茨公式来计算()?x baf 的值的前提是 :()x f 的原函数能够求出。当 ()x f 的原函数不易求出或找不到时 ,希望用一个易于求原函数的函数来近似代替被积函数 ,从而得到定积分的近似计算公式。下文中梯形公式T 就是()?x baf 常用的近似计算公式。 复化梯形公式的主要思想是利用若干小梯形的面积代替原方程的积分，利用微元法，可以求出坐标面上由函数与坐标轴围城的图像的面积的近似值，符合了计算机计算存储的思想。 下面，我们在探讨复化梯形公式的计算规律。 设将求积区间[]b a ,分成n 等份，则一共有1+n 个分点，按梯形公式 ()()[]()()()??????++=+=∑∑-=+-=b f x f a f h x f x f h T n k k k k n k n 1111 0222计算积分值n T ，需要提供1+n 个函数值。 注意，这里n a b h -= 代表步长，分点为kh a x k +=，.,1,0n k ???= 2. 编写用变步长复化辛普森公式计算积分()b a f x dx ? 的程序。 用上面编写的程序计算下列积分并分析计算结果 （1）0cos xdx π ? （2）220cos x x dx （3）?10dx x 程序： function S=bianfuhuasimpson(fx,a,b,eps,M) % 变步长复合simpson 求积公式 % 调用方式： S=fuhuasimpson(@fx,a,b,epsilon) % fx -- 求积函数（函数文件） % a, b -- 求积区间 % eps -- 计算精度 % M--最大允许输出划分数 n=1; h=(b-a)/n; T1=h*(feval(fx,a)-feval(fx,b))/2; Hn=h*feval(fx,(a+b)/2); S1=(T1+2*Hn)/3; n=2*n; % 最好与倒数第三行保持一致（变步长） while n<=M T2=(T1+Hn)/2; Hn=0; h=(b-a)/n; for j=1:n x(j)=a+(j-1/2)*h; y(j)=feval(fx,x(j)); Hn=Hn+y(j); end Hn=h*Hn; S2=(T2+2*Hn)/3; fprintf(' n=%2d S2=%-12.9f S2-S1=%-12.9f\n',n,S2,abs(S2-S1)); if abs(S2-S1) S=S2; % 达到下列条件之一，则运算终止: % (1).abs(S2-S1) 实习七 数值积分的辛普森方法 一、实习目的 1．掌握计算定积分近似值的辛普森方法； 2．理解复化辛普森求积公式。 二、相关知识 抛物线公式（辛普森公式） 将积分区间],[b a 作2n 等分：n i ih a x n a b h i 2,,2,1,0,,2 =+=-=，现在考察由分点22-k x 和k x 2形成的一个小区间],[222k k x x -，（12-k x 为中点），n k ,,2,1 =，在每一个 小区间],[222k k x x -上，作一条抛物线k k k x x y γβα++=2通过三点))(,(2222--k k x f x ， ))(,(1212--k k x f x 和))(,(22k k x f x ，这样就产生关于未知系数k α，k β和k γ的线性方程组 ?????=++=++=++------)() ()(222212122122222222k k k k k k k k k k k k k k k k k k x f x x x f x x x f x x γβαγβαγβα （7-1） 显然上述方程组有唯一解（由高等代数知识知）。 现在，以)(x f y =为顶的曲边梯形用以抛物线k k k x x y γβα++=2为顶的曲边梯形来 代替，其面积 dx x x dx x f k k n k x x k b a k k )()(1 2222γβα++≈∑??=-∑=--=n k k k x x 12226]4)(2)2([){(2222222222222222k k k k k k k k k k k k k x x k x x x x x x γβαγβα++++++++----- }222k k k k k x x γβα+++)]()(4)([6212221222k k k n k k k x f x f x f x x ++-=--=-∑ （7-2） 得抛物线公式，记为n S 2，化简后： {})()(4)(2)(4)(2)(4)(6212432102n n n x f x f x f x f x f x f x f n a b S +++++++-=- 在实际求解数值积分时，我们总是采用成倍加密节点的方法，就抛物线公式而言，若n S 2被认为精度不够，则接着计算n S 4，而精度是否达到要求，又以n n S S 24-是否足够小作为判复化梯形公式及复化辛普森公式的精度比较
利用复化梯形公式、复化simpson 公式计算积分
matlab实现复化梯形公式,复化simpson公式以及romberg积分
龙格库塔积分算法
复化梯形公式和复化Simpson公式word.doc
辛普森求积公式
复合梯形积分和复合Simpson积分计算数值积分
关于辛普森(simpson)公式在线路坐标计算中的应用
数值分析复化Simpson积分公式和复化梯形积分公式计算积分的通用程序
5.2龙格库塔法
辛普森公式
选用复合梯形公式复合Simpson公式计算
复化梯形公式,辛普森公式的matlab程序
验证数值积分求积公式及复合梯形公式程序设计
龙格库塔法例题
复化梯形公式
变步长复化辛普森公式计算积分
数值积分的辛普森方法