数值分析上机实验题参考
数值分析上机实验报告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。
数值分析上机实验题参考
数值分析论文数值积分 一、问题提出选用复合梯形公式，复合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次代数精度。
三、复合求积公式为了提高精度，通常可以把积分区间分成若干子区间（通常是等分），再在每个子区间用低阶求积公式，即复化求积法，比如复化梯形公式与复化辛普森公式。
数值分析第三次上机练习实验报告
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 级数的部分和（有限项）做拟合函数。
(完整word版)数值分析上机作业1-1解析
考虑一个高次的代数多项式
（E1-1）
显然该多项式的全部根为l，2，…，20，共计20个，且每个根都是单重的（也称为简单的）。现考虑该多项式方程的一个扰动
(E1-2)
其中 是一个非常小的数。这相当于是对（E1-1）中 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。
ve=zeros(1,21);
ve(21-Numb)=ess;
root=roots(poly(1:20)+ve);
x0=real(root); y0=imag(root);
plot(x0',y0','*');
disp(['对扰动项 ',num2str(Numb),'加扰动',num2str(ess),'得到的全部根为:']);
ess分别为1e-6,1e-8.1e-10,1e-12的图像如下：
从实验的图形中可以看出，当ess充分小时，方程E.1.1和方程E.1.2的解相差很小，当ess逐渐增大时，方程的解就出现了病态解，这些解都呈现复共轭性质。
（2）将扰动项加到x18上后，ess=1e-009时方程的解都比较准确，没有出现复共轭现象。ess=1e-008时误差与x19(ess=1e-009)时相当，即扰动加到x18上比加到x19小一个数量级。对x8的扰动ess=1000时没有出现复共轭，误差很小；对x的扰动ess=10e10时没有出现复共轭，误差很小。因此，扰动作用到xn上时，n越小，扰动引起的误差越小。
if((Numb>20)|(Numb<0))errordlg('请输入正确的扰动项:[0 20]之间的整数!');return;end
数值分析上机实习题
数值分析上机实习题第2章插值法1. 已知函数在下列各点的值为试⽤四次⽜顿插值多项式）（x p 4及三次样条韩式)(S x (⾃然边界条件)对数据进⾏插值。
⽤图给出(){}10,11,1,0,08.02.0,,x i =+=i x y i i ，）（x p 4及)(x S Python 代码import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.font_manager import FontPropertiesfont_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=12) #求⽜顿n 次均差 def qiujuncha(x,f,n): for i in range(1,n): for j in range(4,i-1,-1):f[j]= (f[j] - f[j-1])/(x[j]-x[j-i]) #根据⽜顿多项式求值 def niudun(x,f,x1): sum = f[0]; tmp = 1;for i in range(1,5): tmp *= (x1-x[i-1]) sum = sum + f[i]*tmp return sum#⽜顿插值画图 def drawPic(x,f):x1 = np.linspace(0.2, 1, 100) plt.plot(x1, niudun(x,f,x1))plt.title(u"⽜顿四次插值",fontproperties=font_set) plt.xlabel(u"x 轴",fontproperties=font_set) plt.ylabel(u"y 轴", fontproperties=font_set) plt.show() def qiu_h(x,h): n = len(x) -1 for i in range(n): print(i)h[i] = x[i+1]-x[i]#⾃然边界条件下的三次样条插值求Mdef qiu_m(h,f,o,u,d):n = len(h)o[0] = 0u[n] = 0d[n] = d[0] = 0a = []for i in range(1,n):u[i] = h[i-1]/(h[i-1]+h[i])for i in range(1,n):o[i] = h[i]/(h[i-1]+h[i])for i in range(1,n-1):d[i] = 6*(f[i+1]-f[i])/(h[i-1]+h[i])t = [0 for i in range(5)]t[0] =2t[1] = o[0]a.append(t)for i in range(1,n):t = [0 for i in range(5)]t[i - 1] = u [i + 1]t[i] = 2t[i + 1] = o [i + 1]a.append(t)t = [0 for i in range(5)]t[n - 1] = u[n]t[n] = 2a.append(t)tmp = np.linalg.solve(np.array(a),np.array(d))m = []for i in range(5):m.append(tmp[i])return m#根据三次条插值函数求值def yangtiao(x1,m,x,y,h,j):returnm[j]*(x[j+1]-x1)**3/(6*h[j])+m[j+1]*(x1-x[j])**3/(6*h[j])+(y[j]-m[j]*h[j]**2/6)*(x[j+1]-x1)/h[j] +(y[j+1]-m[j+1]*h[j]**2/6)*(x1-x[j])/h[j] def main():x = [0.2, 0.4, 0.6, 0.8, 1.0]y = [0.98, 0.92, 0.81, 0.64, 0.38]f = y[:]f1 = y[:]h = [0.2,0.2,0.2,0.2]u = [0 for n in range(5)]d = [0 for n in range(5)]o = [0 for n in range(5)] qiujuncha(x,f,4) qiujuncha(x,f1,2)m = qiu_m(h,f1,o,u,d) x1 = np.linspace(0.2, 0.4, 10)p1= plt.plot(x1, yangtiao(x1,m,x,y,h,0),color='red') x1 = np.linspace(0.4, 0.6, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 1),color='red') x1 = np.linspace(0.6, 0.8, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 2),color='red') x1 = np.linspace(0.8, 1.0, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 3),color='red') x1 = np.linspace(0.2, 1.0, 40)p2 = plt.plot(x1,niudun(x,f,x1),color='green') plt.xlabel(u"x 轴", fontproperties=font_set) plt.ylabel(u"y 轴",fontproperties=font_set) plt.title("三次样条插值和⽜顿插值")plt.legend(labels=[u'三次样条插值',u'⽜顿插值'],prop=font_set,loc="best") plt.show() main()实验结果运⾏结果可得插值函数图（如图1-1），4次⽜顿插值函数）（x p 4和三次样条插值函数)(x S 如下：)6.0(*)4.0(*)2.0(625.0)4.0(*)2.0(*3.098.0)(4-------=x x x x x x x P 98.0)8.0(*)6.0(*)4.0(*)2.0(*20833.0+-----x x x x]4.0,2.0[),2.0(467.4)4.0(9.4)2.0(167.1)(S 3∈-+-+-=x x x x x]6.0,4.0[),4.0(113.4)6.0(6467.4)4.0(575.1)6.0(167.1)(S 33∈-+-+----=x x x x x x ]8.0,6.0[),6.0(2.3)8.0(113.4)6.0(575.1)(S 3∈-+-+--=x x x x x]0.1,8.0[),8.0(9.1)0.1(2.3)(S ∈-+-=x x x x图1-1三次样条插值和⽜顿插值图2.在区间[-1,1]上分别取n = 10,20⽤两组等距节点对龙格函数做多项式插值三次样条插值，对每个n值画出插值函数及图形。
东南大学 数值分析上机题作业 MATLAB版
东南大学数值分析上机题作业MATLAB版上机作业题报告1.Chapter 11.1题目设S N= Nj=2j2−1，其精确值为11311(-- 。
22N N +1（1）编制按从大到小的顺序S N =（2）编制按从小到大的顺序S N =111，计算S N 的通用程序。
++⋯⋯+22-132-1N 2-1111，计算S N 的通用程++⋯⋯+N 2-1(N -1 2-122-1序。
（3）按两种顺序分别计算S 102, S 104, S 106, 并指出有效位数。
（编制程序时用单精度）（4）通过本次上机题，你明白了什么？1.2程序 1.3运行结果1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。
按从小到大的顺序，有效位数分别为：5，6，6。
可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。
当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。
因此，采取从小到大的顺序累加得到的结果更加精确。
2.Chapter 22.1题目（1）给定初值x 0及容许误差ε，编制牛顿法解方程f(x=0的通用程序。
3（2）给定方程f (x =x-x =0, 易知其有三个根x 1*=-3, x 2*=0, x 3*=○1由牛顿方法的局部收敛性可知存在δ>0, 当x 0∈(-δ, Newton 迭代序列收敛+δ 时，于根x2*。
试确定尽可能大的δ。
○2试取若干初始值，观察当x 0∈(-∞, -1, (-1, -δ, (-δ, +δ, (δ, 1, (1, +∞ 时Newton 序列的收敛性以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？2.2程序2.3运行结果（1）寻找最大的δ值。
算法为：将初值x0在从0开始不断累加搜索精度eps ，带入Newton 迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma 值。
数值分析实验(参考答案)
数值分析上机实验学生姓名： 学号： 教师：实验1：1. 实验项目的性质和任务通过上机实验，使学生对病态问题、线性方程组求解和函数的数值逼近方法有一个初步理解。
2．教学内容和要求 1）对高阶多多项式201()(1)(2)(20)()k p x x x x x k ==---=-∏编程求下面方程的解 19()0p x x ε+=并绘图演示方程的解与扰动量ε的关系。
Matlab 程序：x=1:20;y=zeros(1,20); ve=zeros(1,21); plot(x,y,'o') hold on ; pause; for x=1:5ve(2)=10^(-x);e=roots(poly(1:20)+ve);plot(e,'*') hold on pause; end0510152025303540-20-15-10-5051015202）对2~20n =，生成对应的Hilbert 矩阵，计算矩阵的条件数；通过先确定解获得常向量b 的方法,确定方程组 n H x b =最后，用矩阵分解方法求解方程组，并分析计算结果。
Matlab 程序：for n=2:20; H=hilb(n); ca=cond(H,2) x=(1:n)'; b=H*x; [L,U]=lu(H); y=L\b;x1=U\yplot(x,x,'o',x1,x1,'*') hold on pause; end-500-400-300-200-100100200300400500-500-400-300-200-10001002003004005003）对函数 21()[1,1]125f x x x=∈-+的Chebyshev 点 (21)cos()1,2,...,12(1)kk x k n n π-==++编程进行Lagrange 插值，并分析插值结果。
Matlab 程序：function y=lagrangen(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i);s=0;for k=1:nL=1;for j=1:nif j~=kL=L*(z-x0(j))/(x0(k)-x0(j));endends=s+L*y0(k);endy(i)=s;endy;for n=5:20x=-1:0.01:1; y=1./(1+25*x.^2);plot(x,y,'r')hold onk=n+1:-1:1;x0=cos((2*k-1)*pi./(2*(n+1))),y0=1./(1+25*x0.^2);x=-1:0.01:1; y1=lagrangen(x0,y0,x);plot(x0,y0,'o',x,y1), pausehold off end-1-0.8-0.6-0.4-0.200.20.40.60.81-0.200.20.40.60.811.2-1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81。
数值分析上机实例
《数值分析实验报告》实验题目：曲线拟合的最小二乘法一、实验目的：1、掌握拟合曲线的最小二乘法;2、探求拟合函数的选择和拟合精度间的关系。
二、实验要求：1、用最小二乘法进行三次多项式的拟合；2、计算y j与y(t j)的误差，j=1,2…11；3、另外选取一个拟合函数，进行拟合效果比较；4、* 绘制出曲线拟合图。
三、实验问题：在某冶炼过程中，通过实验检测得到的含碳量与时间关系数据如下，试求含碳量y与时间t内在关系的拟合曲线。
四、实验过程与结果：先画出实验数据图形，如下图：应实验要求取拟合函数为三次函数：即：设Y=X(1)+X(2)*t+X(3)*t.^2+X(4)*t.^3;取权函数均为1令x=（0，5,10,15,20,25,30,35,40,45,50）’y=（0，1.27，2.16，2.86，3.44，3.87，4.15，4.37，4.51 ，4.58，4.02）’;x1=（1,1,1,1,1,1,1,1,1,1,1）’;x2=x;x3=x.^2;x4=x.^3;B=[(x1*x1'),(x1*x2'),(x1*x3'),(x1*x4');(x2*x1'),(x2*x2'),(x2*x3'),(x2*x4');(x3*x1'),(x3*x2'),(x3*x3'),(x3*x4');(x4*x1'),(x4*x2'),(x4*x3'),(x4*x4')];C=((x1*y'),(x2*y'),(x3*y'),(x4*y'))’;其矩阵形式是 B*X=C’X=（X（1）, X（2）,X（3）,X（4））’解出得X =（0.0887412587412430.233658896658901-0.0033839160839160.000006790986791）所以得到三次拟合曲线为：Y=0.088741258741243+0.233658896658901*t-0.003383916083916*t.^2+ 0.000006790986791*t.^3 得到曲线拟合图为：求y j与y(t j)均方误差：sum=0sum=(0.088741258741243*x1(j)+0.233658896658901*x2(j)-0.00338 3916083916*x3(j)+ 0.000006790986791*x4(j)-y(j))^2+sumsqrt(sum) j=1,2,3,4,5,6,7,8,9,10,11求得均方误差结果为：ans = 0.317886926283700另外同理选取一个二次多项式进行拟合：得到二次拟合曲线系数为：Q =（ 0.1193006993006950.223947785547786-0.002874592074592）二次拟合曲线为：Y1=0.119300699300695+0.223947785547786*t1-0.002874592074592*t1.^2;二次拟合均方误差为ans =0.324813161284437得到二次拟合图：五、结果分析与讨论将二次多项式和三次多项式的曲线拟合图放在一块比较如下图：可以看出三次拟合要比二次拟合要好一点但并不是很明显，其均方误差比较接近相差-0.006926235000737即二次和三次拟合精度相差较小。
数值分析第五版上机实验答案实验一~实验六
实验一Lagrange插值算法实验目的：掌握拉格朗日（Lagrange）插值算法的基本原理，理解插值基函数的性质，掌握基本的误差概念。
学习用计算机语言编写程序实现算法。
[参考程序]#include "stdio.h"//定义插值节点及所求点数据，根据题目不同而修改double x[]={0.32,0.34,0.36};double y[]={0.314567,0.333487,0.352274};double xx=0.3367;// Lagrange插值算法函数，利用循环计算具有对称性的基函数和最终结果double Lagrange(double xxx,int n){int i;double result=0,temp;for(i=0;i<=n;i++){temp=1;for(int j=0;j<= n;j++){if(j!=i){temp=temp*(xxx-x[j])/(x[i]-x[j]);}}result=result+temp*y[i];}return result;}void main(){int n;printf("Please input n:");scanf("%d",&n);printf("Sin(%f) = %f \n",xx,Lagrange(xx,n));}实验二Newton均差插值算法实验目的：掌握Newton均差插值算法的基本原理，理解均差的概念，掌握均差表的计算方法。
学习用计算机语言编写程序实现算法。
[参考程序]#include "stdio.h"#define N 10double f[N][N];//定义插值节点及所求点数据，根据题目不同而修改double x[]={0.4,0.55,0.65,0.80,0.90,1.05};double y[]={0.41075,0.57815,0.69675,0.88811,1.02652,1.25382};double fx(int i,int j);double S(int start,int end,double xx);main(){int loopi,loopj,n;double result,xx;scanf("%d",&n);scanf("%lf",&xx);for(loopi=0;loopi<=n;loopi++){//零阶均差作为均差表二维数组的第0列f[loopi][0]=y[loopi];}//循环计算均差表中的所有数据for(loopi=1;loopi<=n;loopi++){for(loopj=1;loopj<=loopi;loopj++){f[loopi][loopj]=fx(loopi,loopj);}}result=S(0,n,xx);printf("Result is: %f",result);return 1;}//求均差的函数double fx(int i,int j){if(j==0){return f[i][j];}else{//这种表示方法需要注意两个x的下标return (fx(i,j-1)-fx(i-1,j-1))/(x[i]-x[i-j]);}}//用秦九韶算法计算插值多项式结果double S(int start,int end,double xx){if(start==end){return f[end][end];}else{return (S(start+1,end,xx)*(xx-x[start])+f[start][start]);}}实验三Newton差分插值算法实验目的：掌握Newton差分插值算法的基本原理，理解差分的概念，掌握差分表的计算方法。
数值分析上机实验报告_线性方程组部分实验题1
线性方程组部分实验题11. 问题描述2. 问题分析2.1. Hilbert 矩阵Hilbert 矩阵是一种著名的“坏条件”矩阵，及病态矩阵，该矩阵任意一个元素的微小变化都能引起巨大的误差（如求逆等）。
该矩阵的元素的数学表达式是a(i ，j)=1/(i+j-1)。
这里我们用matlab 自带的产生Hilbert 矩阵的函数hilb(n)计算（n=1、 2......15）。
同时可以用MATLAB 计算。
2.2. Gauss 消去法（LU 分解）和Cholesky 分解方法高斯消去法主要有顺序消去和回代过程，是线性方程最基本的解法之一。
对n n A R ⨯∈，若0detA ≠，0(1,2,,)i i n ∆≠=⋅⋅⋅，则存在唯一的单位下三角形矩阵L 和非奇异上三角形矩阵U 使得A=LU 。
Cholesky 分解方法，也称平方根法。
对n n A R ⨯∈，假设A 对称正定，则存在唯一对角元素为正数的下三角矩阵L ，使得T A=LL 。
2.3. 条件数计算根据条件数定义：1222()n n n cond H H H -= ；12221n i i X x =⎛⎫== ⎪ ⎪⎝⎭∑ 3. 求解方案设计3.1. 线性方程组的求解及误差在这里为了比较，选择n=2、5、8、10、15 五种维数进行计算分析。
按照题意，我们先设x 为分量全1的向量，求出b ，然后将H 和b 作为已知量，求x ，与设定的参考解对比。
Gauss 消去法、Cholesky 分解法的Matlab 实现程序如下：----------------------------------------------------------------------------------------------------------%Guass 法***************************************************************************①Guass.m function guass(n)n=input('请输入系数矩阵的维数n='); %输入nH=hilb(n); %构造Hilbert 矩阵 x_exact=ones(n ，1); %构造x 列向量 b=H*x_exact; %构造b 向量 x=Doolittle(H ，b)a=input('是否继续，继续请按1，结束请按0:'); if a==1clear;clc;close all; guass(n); else end②Doolittle.mfunction x=Doolittle(A ，b) % LU 分解求Ax=b 的解 N=size(A); n=N(1);L=eye(n ，n);%L 的对角元素为1 U=zeros(n ，n); %U 为零矩阵 U(1，1:n)=A(1，1:n)%U 第一行L(1:n，1)=A(1:n，1)/U(1，1)%L第一列for k=2:nfor i=k:nU(k，i)=A(k，i)-L(k，1:(k-1))*U(1:(k-1)，i)endfor j=(k+1):nL(j，k)=(A(j，k)-L(j，1:(k-1))*U(1:(k-1)，k))/U(k，k)endendy=solveDownTriangle(L，b);%调用下三角系数矩阵求解函数x=solveUpTriangle(U，y);%调用上三角系数矩阵求解函数③solveDownTriangle.mfunction x=solveDownTriangle(A，b)%求下三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=1:nif (i>1)s=A(i，1:(i-1))*x(1:(i-1)，1);elses=0;endx(i，1)=(b(i)-s)/A(i，i);end④solveUpTriangle.mfunction x=solveUpTriangle(A，b)% 求上三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=n:-1:1if (i<n)s=A(i，(i+1):n)*x((i+1):n，1);elses=0;endx(i，1)=(b(i)-s)/A(i，i);end%Cholosky分解方法******************************************************************** ---------------------------------------------------------------------------------------------------------- %利用对称正定矩阵之Cholesky分解求解线性方程组Ax=bfunction x=Chol_Solve(A，b)clear;clc% A=input('输入A=');% b=input('输入b=');n=input('请输入系数矩阵的维数n=');%定义hilbert矩阵A=hilb(n);x_exact=ones(n，1);b=A*x_exact;n=length(b);l=Cholesky(A);x=ones(n，1);y=ones(n，1);y(1)=b(1)/l(1，1);%求Ly=bfor i=2:nz=0;for k=1:i-1z=z+l(i，k)*y(k);endy(i)=(b(i)-z)/l(i，i);endx(n)=y(n)/l(n，n);%求L'x=yfor i=n-1:-1:1z=0;for k=i+1:nz=z+l(k，i)*x(k);endx(i)=(y(i)-z)/l(i，i);endxwc=x-x_exact %误差计算a=input('是否继续，继续请按1，结束请按0:');if a==1Chol_Solve(A，b);else end%对称正定矩阵之Cholesky分解法：function L=Cholesky(A)% n=input('请输入系数矩阵的维数n=');%定义hilbert矩阵% A=hilb(n);n=length(A);L=zeros(n);for k=1:ndelta=A(k，k);for j=1:k-1delta=delta-L(k，j)^2;endif delta<1e-10return;endL(k，k)=sqrt(delta);for i=k+1:nL(i，k)=A(i，k);for j=1:k-1L(i，k)=L(i，k)-L(i，j)*L(k，j);endL(i，k)=L(i，k)/L(k，k);endend3.2.条件数的计算根据2.3叙述原理，写出条件数的计算程序如下：%tiaojianshu.mm=input ('input m:=') ;N=[1:m];for i=1:mn=N(i);H=hilb(n);k=cond(H);disp('矩阵的阶数')disp(n)disp('矩阵')disp(H)disp('矩阵的条件数')disp(k)end4.结果与误差分析4.1.线性方程组计算结果分别用Gauss法和Cholesky分解法计算线性方程组的结果，并进行误差分析结果如下：4.2.条件数计算结果4.3. 结果分析取不同的n 值，得到如下结果:对于Guass 法，可以看出来，随着n 的增大，求解结果误差变大，这是因为随着n 增大，系数矩阵的条件数变大，微小的扰动就容易造成很大的误差。
- 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这只是其中程序的一部分，从结果可以看出，是比较接近精确值的。