数值分析上机实验题参考
数值分析上机实验报告3
实验报告三题目：函数逼近——曲线拟合目的：掌握曲线拟合基本使用方法数学原理：[P,S]=polyfit(x,y,3)其中x，y为取样值，3为得出的结果的最高次数。
P为对应次数的系数，S为误差值向量，其中x,y是等长的向量，P是一个长度为m+1的向量。
结果分析和讨论：23.观察物体的运动，得出时间t与距离s的关系如表，求运动方程。
t=[0,0.9,1.9,3.0,3.9,5.0];s=[0,10,30,50,80,110];[P,S]=polyfit(t,s,5)P =-0.5432 6.4647 -26.5609 46.1436 -13.2601 -0.0000S =R: [6x6 double]df: 0normr: 1.2579e-012所以得到方程为：-13.2601x46.1436x-26.5609x6.4647x-0.5432x2345++24.在某化学反应堆里，根据实验所得分解物的质量分数y与时间t的关系，用最小拟合求y=F（t）；>> x=0:5:55;y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.62,4.64];>> [P,S]=polyfit(x,y,5)P =0.0000 -0.0000 0.0002 -0.0084 0.2851 0.0082S =R: [6x6 double]df: 6normr: 0.0487所以得到方程为：0082.02851.00084.00002.023++-xxx结论：在23题中计算的结果误差为4.5769，而在24中计算的结果误差为0.0487，说明对于曲线拟合来说，总会有误差，因为取样点并不是都过拟合的曲线的。
数值分析上机实践题3-2013
数值分析上机实践题第三次上机题目(二分法)第一组:组长：李龙宇，组员：杜彦霖，胡朋，黄湘云，雷盛华，李伟元 用二分法求方程010423=-+x x , ]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第二组: 组长：王宇彬，组员：马泽川，权涛涛，师楠？？，路世伦，仲晓磊 用二分法求方程04442234=++--x x x x ,]0,2[-∈x 的近似根,要求根精确到 510- ,并求二分次数.第三组: 组长： 薛原 ，组员：谢胜权，杨帆，王正奇，肖特，张锡云 用二分法求方程04442234=++--x x x x ,]2,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第四组: 组长：柴春晓 ，组员： 韩静兰，李金慧，刘从，马超群，孟凯悦 用二分法求方程04442234=++--x x x x ,]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第五组: 组长：龙纯鹏，组员：代喜，白鑫，鲍亚强，周邦安，张佳伟 用二分法求方程02=--x x ,]1,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第六组: 组长：何关瑶 ，组员：纪伟亮，侯佳意，李济言，李振华，马文磊 用二分法求方程06cos 2=-++-x e e x x ,]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第七组: 组长：杨钦 ，组员： 王凌宇，吴凯杰，薛小龙，袁权炜，师俊峰 用二分法求方程0232=-+-x x e x ,]1,0[∈x 的近似根,要求根精确到 510- ,并求二分次数.第八组: 组长：汪芳 ，组员：张学利，周幸茹，李雨珏，张飞用二分法求方程05.0cos 2=++x x π,]5.1,5.0[∈x 的近似根, 要求根精确到510- ,并求二分次数.第九组: 组长：刘永鸿 ，组员：黄尚政，李超，郭新磊，何奎奎用二分法求 15 的近似根,要求根精确到 510- ,并求二分次数.第十组: 组长：杨吉望 ，组员：龙力，任金雄，王亮，王文强，谢丁波 用二分法求方程 325 的近似根,要求根精确到 510- ,并求二分次数.第十一组: 组长： 张国强，组员： 赵奇，袁硕，郭凯旋，于沛生，鲍宏雷 用二分法求方程 ,05.0c o s 2=++x x π在]2,0[内的近似根,要求根精确到 510- ,并求二分次数.第十二组: 组长：苏映雪 ，组员： 邓晓庆，钟桂平，崔楚轩，高鹏程 用二分法求方程 0797*******=-+--x x x x 的靠近x=2的近似根,要求根精确到 510- ,并求二分次数.备用题:第一组:用二分法求方程 016=--x x , ]2,1[∈x 的近似根,要求根精确到 510- ,并求二分次数.第二组:用二分法求方程 0t a n =-x x ,]5.4,4[∈x 的近似根,要求根精确到 510- ,并求二分次数.补充知识MATLAB 中自带的求根函数:1. roots :求解多项式P(x)=0的根可以用此语句, 输入多项式P(x)的系数(按降幂排列), 输出为P(x)=0的全部根;例如:要求013178)(39=+-+=x x x x P 的根,可以用以下语句:>> fa =[8,0,0,0,0,0,17,0,-3,1]>> gen= roots(fa)运行后输出全部根.2. fsolve: 求解超越方程f(x)=0的根可以用此语句(也可以解多项式方程,但计算量较大), 输入多项式P(x)的系数(按降幂排列), 输出为P(x)=0的全部根调用格式: X = fsolve(F,X0)其中输入函数F(x)的M文件名和解X的初始值X0，X0可以是矩阵或向量。
数值分析上机实习题
2019-2020 第1学期数值分析上机实习题总目标：会算，要有优化意识。
（以下程序要求以附件1例题代码格式给出）1. 对给定的线性方程组Ax b =进行迭代求解。
（1）给出Jacobi 迭代的通用程序。
（2）给出Gauss-Seidel 迭代的通用程序。
调用条件：系数矩阵A ，右端项b ，初值0x ，精度要求ε。
输出结果：方程组的近似解。
给定线性方程组211122241125x --⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪---⎝⎭⎝⎭，和122711122215x -⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，取初值0x 为0， 分别利用Jacobi 迭代和G-S 迭代进行求解，观察并解释其中的数学现象。
2. 利用紧凑格式（即直接分解法或逐框运算法）对给定的矩阵A 进行Doolittle 分解，并用其求线性方程组的解。
调用条件：矩阵A 。
输出结果：单位下三角矩阵L 和上三角矩阵U 。
给定矩阵1112A ⎛⎫= ⎪⎝⎭，利用以下算法：1）将A 作Doolittle 分解11A LU =，2）令211A U L =，并对2A 作Doolittle 分解222A L U =，3）重复2）的过程令11n n n A U L --=，并对n A 作Doolittle 分解n n n A L U =，2,3,4,n =， 观察n L ，n U ，n A 的变化趋势，思考其中的数学现象。
3. 给定函数21(),12511f x x x -≤+≤=，取164,8,n =，用等距节点21,i i n x =-+ 0,1,,1i n =+对原函数进行多项式插值和五次多项式拟合，试画出插值和拟合曲线，并给出数学解释。
4. 给出迭代法求非线性方程()0f x =的根的程序。
调用条件：迭代函数()x ϕ，初值0x输出结果：根的近似值k x 和迭代次数k给定方程32()10f x x x =--=，用迭代格式1k x +=0 1.5x =附近的根，要使计算结果具有四位有效数字，利用估计式*1||1||k k k L x x x x L -≤---，或估计式*10||1||kk L x x x x L-≤--来判断需要的迭代次数，分别需要迭代多少次？两者是否有冲突？5. 利用数值求积算法计算()ba f x dx ⎰。
数值分析上机题目
数值分析上机题目4(总21页) --本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--实验一实验项目：共轭梯度法求解对称正定的线性方程组 实验内容：用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k x x --∞-<时停止计算。
编制程序：储存m 文件function [x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;while rho>10^(-11) & k<1000 k=k+1; if k==1 p=r; elsebeta=rho/rho1; p=r+beta*p; end w=A*p;alpha=rho/(p'*w); x=x+alpha*p; r=r-alpha*w; rho1=rho;rho=r'*r; end运行程序： clear,clcA=[2 -1 0 0;-1 3 -1 0;0 -1 4 -1;0 0 -1 5]; b=[3 -2 1 5]'; [x,k]=CGmethod(A,b)运行结果： x =(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵， A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。
编制程序：储存m 文件function [x,k]=CGmethod_1(A,b) n=length(A);x(1:n,1)=0;r=b-A*x;r1=r; k=0;while norm(r1,1)>10^(-7)&k<10^4 k=k+1; if k==1 p=r; elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p; end r=r1; w=A*p;alpha=(r'*r)/(p'*w); x=x+alpha*p; r1=r-alpha*w; end运行程序： clear,clc n=1000; A=hilb(n); b=sum(A')';[x,k]=CGmethod(A,b)实验二1、 实验目的：用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。
东南大学数值分析上机题答案
东南⼤学数值分析上机题答案数值分析上机题第⼀章17.（上机题）舍⼊误差与有效数设∑=-=Nj N j S 2211，其精确值为）111-23（21+-N N 。
（1）编制按从⼤到⼩的顺序1-1···1-311-21222N S N +++=，计算N S 的通⽤程序；（2）编制按从⼩到⼤的顺序121···1)1(111222-++--+-=N N S N ，计算NS 的通⽤程序；（3）按两种顺序分别计算210S ,410S ,610S ,并指出有效位数（编制程序时⽤单精度）；（4）通过本上机题，你明⽩了什么？解：程序：（1）从⼤到⼩的顺序计算1-1···1-311-21222N S N +++=：function sn1=fromlarge(n) %从⼤到⼩计算sn1format long ; sn1=single(0); for m=2:1:nsn1=sn1+1/(m^2-1); end end（2）从⼩到⼤计算121···1)1(111222-++--+-=N N S N function sn2=fromsmall(n) %从⼩到⼤计算sn2format long ; sn2=single(0); for m=n:-1:2sn2=sn2+1/(m^2-1); end end（3）总的编程程序为： function p203()clear allformat long;n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为sn fprintf('精确值为%f\n',sn);sn1=fromlarge(n);fprintf('从⼤到⼩计算的值为%f\n',sn1);sn2=fromsmall(n);fprintf('从⼩到⼤计算的值为%f\n',sn2);function sn1=fromlarge(n) %从⼤到⼩计算sn1 format long;sn1=single(0);for m=2:1:nsn1=sn1+1/(m^2-1);endendfunction sn2=fromsmall(n) %从⼩到⼤计算sn2 format long;sn2=single(0);for m=n:-1:2sn2=sn2+1/(m^2-1);endendend运⾏结果：从⽽可以得到N值真值顺序值有效位数2 100.740050 从⼤到⼩0.740049 5从⼩到⼤0.740050 64 100.749900 从⼤到⼩0.749852 3从⼩到⼤0.749900 66 100.749999 从⼤到⼩0.749852 3从⼩到⼤0.749999 6（4）感想：通过本上机题，我明⽩了，从⼩到⼤计算数值的精确位数⽐较⾼⽽且与真值较为接近，⽽从⼤到⼩计算数值的精确位数⽐较低。
数值分析上机实习题及答案.docx
方詡文金兴:爭［数值分析］2017-2018第二学期上机实习题1：编程计算亍丄，其中C= 4. 4942x10307,给出并观察计算结心C"果，若有问题，分析之。
解:mat lab 编程如下:E) funct ion diy i ti formatlong g;n 二input C 输入ii 值= c= 4.4942E307; sum 二0; s 二 0;E3 for i = l:n s = l/ (c*i);>> diyiti 输入n 值:10 104.6356e-308 >> diyiti输入ri 值:1001004.6356e-308 >> diyiti 输入n 值:1000 10004.6356e-308 >> diyiti揄入n 值* 1000001000004・ 6356e-308 >> diyiti输入n 值；1000000001000000004.6356e-308图二：运行结果Mat lab 中，forma t long g 对双精度，显示15位定点或浮点格式，由上图 可知，当输入较小的n 值5分别取10, 100, 1000, 100000, 100000000）的时候, 结果后面的指数中总是含有e-308,这和题目中的C 值很相似，我认为是由于分 母中的C 值相对于n 值过大，出现了 “大数吃小数”的现彖，这是不符合算法原 则的。
2：利用牛顿法求方程-1^ = 2于区间241的根，考虑不同初值下牛顿法的收敛情况。
解：牛顿法公式为:利用mat lab 编程function di2ti21 3i=l ；2 2.85208156699784 xO 二input （'输入初值x0：‘ ）;A 二[i x0];3 2.55030468822809 t=x0+ (x0-log (xO) -2) /(1-1/xO) ; %迭代函数4 1.91547247100476 三 while (abs (t _x0)>0.01)i=i+l; 5 0.37867158538991 xO 二 t; 6 0.774964959780275 A = [A;i xO];t =x0+(x0-log(xO)-2)/(1-1/xO): 7 4.11574081641933 cnd| 8 5.04162436446126 disp (A);96.81782826645596当输入初值二3的时候并不能收敛。
数值分析上机题3
数值分析上机题目3实验一1．根据Matlab 语言特点，描述Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法。
2．编写Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
3．给定2020⨯∈R A 为五对角矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------321412132141412132141412132141412132141213O O O O (1)选取不同的初始向量)0(x 及右端面项向量b ，给定迭代误差要求，分别用编写的Jacobi 迭代法和Gauss-Seidel 迭代法程序求解，观察得到的序列是否收敛？若收敛，通过迭代次数分析计算结果并得出你的结论。
(2)用编写的SOR 迭代法程序， 对于(1)所选取的初始向量)0(x 及右端面项向量b 进行求解，松驰系数ω取1<ω<2的不同值，在5)1()(10-+≤-k k x x 时停止迭代，通过迭代次数分析计算结果并得出你的结论。
实验11、 根据MATLAB 语言特点，描述Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法。
2、 编写Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
Jacobi 迭代法function [x1,k]=GS_2(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=D\((L+U)*x0+b);endk=kx=x1Gauss-Seidel迭代法function [x1,k]=GS_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=(D-L)\U*x0-D\b;endk=kx=x1SOR迭代法function [x1,k]=SOR_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0;w=0.96;while norm(x1-x0,1)>10^(-7)&k<100k=k+1;x0=x1;x1=(D-w*U)\(((1-w)*D+w*L)*x0+w*b);endk=kx=x13、采用Jacobi迭代法，Gauss-Seidel迭代法求解五对角矩阵clear,clcA=diag(3*ones(20,1))+diag((-0.5)*ones(19,1),-1)+diag((-0.5)*ones(19,1 ),1)+diag((-0.25)*ones(18,1),-2)+diag((-0.25)*ones(18,1),2);b=sum(A')';[x1,k1]=Jacob_h(A,b)[x2,k2]=GS_h(A,b)运行结果：两种方法都收敛，k1=27,k2=13。
数值分析上机实验
刘力辉 2010210804011.“画圆为方”问题也是古希腊人所提出几何三大难题中的另一个问题。
即求作一个正方形，使其面积等于已知圆的面积。
不妨设已知圆的半径为 R = 1，试用数值试验显示“画圆为方”问题计算过程中的误差。
（1）MATLAB 程序：y=pi^(1/2); % to generate 15-bit value of square root of pi b=1; d=1; for k=1:8 b=b*10;d=d/10; % b and d combined to control the digit of x x=d*fix(b*y); s(k)=x^3; l(k)=x; endformat long [l',s']（2）误差分析：2． 设，I x xd x n n=+⎰501由 x n = x n + 5 x n – 1 – 5 x n – 1 可得递推式 I n = – 5I n – 1 + 1/ n（1）从 I 0 尽可能精确的近似值出发，利用递推公式：I I nn n =-+-511 ( n = 1,2,…20)计算从 I 1 到 I 20 的近似值；（2）从 I 30 较粗糙的估计值出发，用递推公式：I I nn n -=-+11515 ( n= 30,29, …, 3, 2 )计算从I 1 到 I 20 的近似值；(3) 分析所得结果的可靠性以及出现这种现象的原因。
I 0 =dx x⎰+1051=ln (5+x )10|=ln 6-ln 5 所以I0≈0.18232155679395format longI0=log2(6)/log2(exp(1))-log2(5)/log2(exp(1)) % calculate the value of I0=ln6-ln5 for n=1:20I0=-5*I0+1/n; % recycling equation between I(n+1) and I(n) s(n)=I0; end s'则计算结果为：表1I1 0.0883922160302300 I11 0.0140713362538500 I2 0.0580389198488700 I12 0.0129766520640700 I3 0.0431387340890000 I13 0.0120398166027400 I4 0.0343063295550100 I14 0.0112294884148600 I5 0.0284683522249700 I15 0.0105192245923700 I6 0.0243249055418100 I16 0.0099038770381400 I7 0.0212326151481100 I17 0.0093041442210800 I8 0.0188369242594600 I18 0.0090348344501700 I9 0.0169264898137900 I19 0.0074574066965100 I10 0.0153675509310500I20 0.0127129665174600从计算的数据看出I 20=0.0127129665174600 > I 19=0.0074574066965100又I n 的积分范围为0～1，所以应该有I n >I n+1。
数值分析第三次上机练习实验报告
2 1 , k 0,1, 2...n. 这里唯一需要注意的就是 n 1 Pn ( xk ) Pn'1 ( xk )
还需要一个多项式微分运算（使用的是polyder函数） ，求解出这两项之后，剩下的 就是 I n ( f ) 2.
A
k 0
n
k
f ( xk ) 。
Romberg积分公式计算 函数主体是Romberg.m文件,以下只分析这个函数文件的结构: 第一部分是赋初值,对于积分表T还赋值T(1,1):
实现方法说明
程序语言采用Matlab语言，运行环境为Matlab R2009b。
注：主程序文件是First.m和Second.m，其他都是对应的自编函数文件。
1.
复化Guass-Legendre求积公式运算 这个程序主体是First.m文件，在里面定义了变量x和函数f f=sym(100*sin(10*(2+x)^(-1))*(2+x).^(-2));由于Legendre多项式的使用条件是[-1,1]， 因此对积分区间做了变换，函数体本身也就变成了f（x）=(100*sin(10/(x + 2)))/(x + 2)^2。其他都是调用自己编写的函数计算。 ① Legendre多项式生成函数：Legendre.m
2/5
水利系
2008010249
程国安
这个函数主要用于生成 Legendre 多项式，比较简短，最重要的就是一个迭代 关系 P=((2*n-1)*x*Legendre(n-1)-(n-1)*Legendre(n-2))/(n);最后输出一个多项式， 如果需要整合，另外加一句 expand（P）即可。 ② Guass-Legendre积分函数：Guass-Legendre.m 这个函数的相关输入输出说明见函数文件注释。 积分实现主要需要求出 f ( xk ) 和 Ak ， 对于 f ( xk ) ， 重点就是建立内联函数 f ( x) 和 变量序列 xk ，变量序列 xk 主要使用利用Legendre求出对应次数的Legendre多项式 的零点（使用sym2poly求出系数，之后利用roots求零点） 。 对于 Ak ，利用 Ak
数值分析上机实验6
数值分析上机实验6根据表中数据，预测公元2000年时的世界人口。
问题分析与数学模型设人口总数为N(t)，根据人口理论的马尔萨斯模型，采用指数函数N(t) = e a + b t=+，令对数据进行拟合。
为了计算方便，将上式两边同取对数，得ln N a bty = ln N或N = e y变换后的拟合函数为y(t) = a + b t根据表中数据及等式 k k( 1，2，……，9)可列出关于两个未知数、b的9个方程的超定方程组（方程数多于未知数个数的方程组）a + t j b = y j(j= 1，2， (9)可用最小二乘法求解。
算法与数学模型求解算法如下：第一步：输入人口数据，并计算所有人口数据的对数值；第二步：建立超定方程组的系数矩阵，并计算对应的正规方程组的系数矩阵和右端向量；第三步：求解超定方程组并输出结果：a，b；第四步：利用数据结果构造指数函数计算2000年人口近似值N(2000)，结束。
MATLAB程序t=1960:1968;t0=2000;N=[29.72 30.61 31.51 32.13 32.34 32.85 33.56 34.20 34.83];y=log(N);A=[ ones(9,1), t' ];d=A\ y' ;a=d(1),b=d(2)N0=exp(a+b*t0)x=1960:2001;yy=exp(a+b*x);plot(x,yy,t,N,'o',2000,N0,'o')计算结果为a =－3，b =6N (2000)所以取五位有效数，可得人口数据的指数拟合函数t e t N 0186.00383.33)(+-=经计算得2000年人口预测值为： (亿)。
例2．温度数据的三角函数拟合问题 洛杉矶郊区在11月8日的温度记录如下在不长的时期内，气温的变化常以24小时为周期，考虑用Fourier 级数的部分和（有限项）做拟合函数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析论文数值积分 一、问题提出选用复合梯形公式，复合Simpson 公式，Romberg 算法，计算I = dx x ⎰-4102sin 4 ()5343916.1≈II =dx x x ⎰1sin ()9460831.0,1)0(≈=I fI = dx xe x⎰+1024 I =()dx x x ⎰++1211ln 二、要求编制数值积分算法的程序；分别用两种算法计算同一个积分，并比较其结果；分别取不同步长()/ a b h -=n ，试比较计算结果（如n = 10, 20等）； ﹡给定精度要求ε，试用变步长算法，确定最佳步长﹡。
三、目的和意义深刻认识数值积分法的意义； 明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题引言一、数值求积的基本思想实际问题当中常常需要计算积分，有些数值方法。
如微分方程和积分方程的求解，也都和积分计算相联系。
依据人们熟悉的微积分基本原理，对于积分I=⎰a b f(x)dx,只要找到被积函数f(x)和原函数F(x),便有下列牛顿-莱布尼茨公式：I=⎰a b f(x)dx=F(b)-F(a).但实际使用这种求积方法往往有困难，因为大量的被积函数，诸如x xsin，2xe-等，其原函数不能用初等函数表达，故不能用上述公式计算。
另外，当f(x)是由测量或数值计算给出的一张数据表时，牛顿-莱布尼茨公式也不能直接运用，因此有必要研究积分的数值计算问题。
二、数值积分代数精度数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对“尽可能多”的函数准确成立，就提出了所谓代数精度的概念。
如果某个求积公式对次数不超过m的多项式均能准确成立，但对m+1次多项式就不能准确成立，则称该求积公式具有m次代数精度。
三、复合求积公式为了提高精度，通常可以把积分区间分成若干子区间（通常是等分），再在每个子区间用低阶求积公式，即复化求积法，比如复化梯形公式与复化辛普森公式。
在区间不大时,用梯形公式、辛普森公式计算定积分是简单实用的,但当区间较大时,用梯形公式、辛普森公式计算定积分达不到精确度要求. 为了提高计算的精确度，我们将[a,b] 区间n等分，在每个小区间上应用梯形公式、辛普森公式计算定积分, 然后将其结果相加，这样就得到了复化梯形公式和复化辛普森公式。
1. 复化梯形公式及余项将积分区间[a,b]n 等分,分点k x a kh =+，b ah n -=, 0,1,.....,,k n = 对每个子区间上采用梯形公式,则得1111n 0()()[()()]2k kn n bx k k ax k k h I f x dx f x dx f x f x R +--+=====++∑∑⎰⎰（f ）记11101[()()][()2()()]22n n n k k k k k h hT f x f x f a f x f b --+===+=++∑∑称为复化梯形公式. 其余项为又因为在区间[a,b]上连续，由连续函数的性质知，在区间[a,b]上存在一点η，于是称为复化梯形公式的余项。
2．复化辛普森公式及余项将区间[a,b]分为n 等份在,每个子区间,1[]k k x x +应用辛普森公式，若记1212k k x x h+=+ 则得 [])()()(4)(6)()(112/111f Rn x f x f x f h dx x f dx x f a bI n k k k k x x n k k k+++===∑⎰∑⎰-=++-=+记[]⎥⎦⎤⎢⎣⎡+++=++=∑∑∑-=+-=+-=++111102/11012/1)()(2)(4)(6)()(4)(6n k k n k k n k k k k b f x f x f a f h x f x f x f h Sn 称为复化辛普森公式。
当f(x)在[a,b]上有连续的四阶导数时，则复化辛普森公式余项推导如下：在区间上辛普森公式的误差已知为所以，在区间[a,b]上公式的误差为又因为在[a,b]上连续, 由连续函数的性质知，在区间[a,b]上存在一点η，于是称为复化辛普森公式的余项。
(二)函数程序编写及结果显示分析在本次函数程序编写中，区间[a,b]默认最多等分为50次。
在具体运行过程，根据题目要求进行具体等分，比如10次，20次等等。
（1）I =dxx ⎰-412sin41、函数程序编制A、总程序#include<stdio.h>#include<math.h>main(){float a,b,h,Tn,Sn,I;float x,y,z;int n,k;printf("请输入区间值a，b \n");printf("a=");scanf("%f",&a);printf("b=");scanf("%f",&b);for(n=1;n<=50;n++){h=(b-a)/n; /*定义步长*/printf("h=%f\n",h);Tn=(h/2)*sqrt(4-sin(a)*sin(a))+(h/2)*sqrt(4-sin(b)*sin(b));/*编写复化梯形公式*//*定义初始值*/for(k=1;k<=n-1;k++){x=k*h;y=(k+1)*h;Tn=Tn+2*(h/2)*sqrt(4-sin(x)*sin(x));}printf("n=%d,Tn=%f\n",n,Tn);/*编写复化辛普森公式*//*定义初始值*/Sn=0;for(k=0;k<=n-1;k++){x=k*h;y=x+(1/2)*h;z=(k+1)*h;Sn=Sn+(h/6)*((sqrt(4-(sin(x)*sin(x))))+4*(sqrt(4-(sin(y)*sin(y))))+(sqrt(4-(sin(z)*sin(z)))));}printf("n=%d,Sn=%f\n",n,Sn);}}2、结果显示A、全部结果显示h=0.00713n=35,Tn=0.498711n=35,Sn=0.498747h=0.006944n=36,Tn=0.498711n=36,Sn=0.498746h=0.006757n=37,Tn=0.498711n=37,Sn=0.498745h=0.006579n=38,Tn=0.498711n=38,Sn=0.498744h=0.006410n=39,Tn=0.498711n=39,Sn=0.498744h=0.006250n=40,Tn=0.498711n=40,Sn=0.498743h=0.006098n=41,Tn=0.498711n=41,Sn=0.498742h=0.005952n=42,Tn=0.498711n=42,Sn=0.498741这只是其中程序的一部分，从结果可以看出，是比较接近精确值的。
B、当取n=20和n=45时，复化梯形公式与复化辛普森公式积分结果显示N=20时，Tn和Sn如下所示：a=0b=0.25n=20h=0.12500n=20,Tn=0.498710n=20,Sn=0.498774N=45时，Tn和Sn如下所示：a=0b=0.25n=45h=0.005556n=45,Tn=0.498711n=45,Sn=0.4987393、结果分析由全部的结果显示可以看出，最后的结果是比较接近准确值I，由于试题给出的结果I=1.5343916，经过推导，程序编制，图像分析，题目给出的准确值I=1.5343916是有误的，准确值大概在0.5左右。
当n=20与n=45时，比较上面两个结果Tn和Sn，它们的计算量基本相同，然而精度确差别不是很大，复化梯形法的结果Tn=0.498710与0.498711，而复化辛普森法的结果Sn=0.0.498774与0.498739。
四个数据基本是都具有4个有效数字。
（2）I =dxxx⎰1sin()9460831.0,1)0(≈=If1、函数程序编制A、总程序#include<stdio.h>#include<math.h>main(){float a,b,h,Tn,Sn,I;float x,y,z;int n,k;printf("请输入区间值a，b\n");printf("a=");scanf("%f",&a);printf("b=");scanf("%f",&b);for(n=1;n<=50;n++){h=(b-a)/n; /*定义步长*/printf("h=%f\n",h);/*编写复化梯形公式*//*定义初始值*/Tn=(h/2)*(1+(h/2)*(sin(b)/b));for(k=1;k<=n-1;k++){x=k*h;y=(k+1)*h;Tn=Tn+2*(h/2)*(sin(x)/x);}printf("n=%d,Tn=%f\n",n,Tn);/*编写复化辛普森公式*//*定义初始值*/Sn=0;for(k=1;k<=n;k++){x=k*h;y=x+(1/2)*h;z=(k+1)*h;Sn=Sn+(h/6)*((sin(x)/x)+4*(sin(y)/y)+(sin(z)/z)); }printf("n=%d,Sn=%f\n",n,Sn);}}2、结果显示由于n的次数过小，不够接近准确值，因此将n扩大，结果如下所示：A、全部结果显示h=0.000502n=1993,Tn=0.945872n=1993,Sn=0.946029h=0.000502n=1994,Tn=0.945870n=1994,Sn=0.946029h=0.000501n=1995,Tn=0.945871n=1995,Sn=0.946030h=0.000501n=1996,Tn=0.945872n=1996,Sn=0.946030h=0.000501n=1997,Tn=0.945873n=1997,Sn=0.946030h=0.000501n=1998,Tn=0.945873n=1998,Sn=0.946030h=0.000500n=1999,Tn=0.945873n=1999,Sn=0.946031h=0.000500n=2000,Tn=0.945873n=2000,Sn=0.946030这只是其中程序的一部分，从结果可以看出，是比较接近精确值的。