数值计算方法大作业
数值计算方法练习题
数值计算方法练习题习题一1. 下列各数都是经过四舍五入得到的近似数，试指出它们有几位有效数字以及它们的绝对误差限、相对误差限。
（1）；（2）；（3）；（4）；（5）；（6）；（7）；2. 为使下列各数的近似值的相对误差限不超过，问各近似值分别应取几位有效数字？3. 设均为第1题所给数据，估计下列各近似数的误差限。
（1）；（2）；（3）4. 计算，取，利用下列等价表达式计算，哪一个的结果最好？为什么？（1）；（2）；（3）（4）5. 序列满足递推关系式若（三位有效数字），计算时误差有多大？这个计算过程稳定吗？6. 求方程的两个根，使其至少具有四位有效数字（要求利用。
7. 利用等式变换使下列表达式的计算结果比较精确。
（1）；（2）（3）；（4）8. 设，求证：（1）（2）利用（1）中的公式正向递推计算时误差增大；反向递推时误差函数减小。
9.设x>0,x*的相对误差为δ，求f(x)=ln x的误差限。
10.下列各数都是经过四舍五入得到的近似值，试指出它们有几位有效数字，并给出其误差限与相对误差限。
11.下列公式如何才比较准确？（1）（2）12.近似数x*=0.0310,是位有数数字。
13.计算取，利用式计算误差最小。
四个选项：习题二1. 已知，求的二次值多项式。
2. 令求的一次插值多项式，并估计插值误差。
3. 给出函数的数表，分别用线性插值与二次插值求的近似值，并估计截断误差。
0.4 0.5 0.6 0.7 0.80.38942 0.47943 0.56464 0.64422 0.717364. 设，试利用拉格朗日余项定理写出以为节点的三次插值多项式。
5. 已知，求及的值。
6. 根据如下函数值表求四次牛顿插值多项式，并用其计算和的近似值。
X 1.615 1.634 1.702 1.828 1.921F (x) 2.41450 2.46459 2.65271 3.03035 3.340667. 已知函数的如下函数值表，解答下列问题（1）试列出相应的差分表；（2）分别写出牛顿向前插值公式和牛顿向后插值公式。
数值分析大作业一
数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。
再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。
2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。
3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。
4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。
二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。
(完整版)数值计算方法试题及答案
数值计算方法试题一一、 填空题（每空1分，共17分）1、如果用二分法求方程043=-+x x 在区间]2,1[内的根精确到三位小数，需对分（ ）次。
2、迭代格式)2(21-+=+k k k x x x α局部收敛的充分条件是α取值在（ ）。
3、已知⎪⎩⎪⎨⎧≤≤+-+-+-≤≤=31)1()1()1(2110)(233x c x b x a x x x x S 是三次样条函数，则a =( )，b =（ ），c =（ ）。
4、)(,),(),(10x l x l x l n 是以整数点n x x x ,,,10 为节点的Lagrange 插值基函数，则∑==nk kx l0)(( )，∑==nk k jk x lx 0)(( )，当2≥n 时=++∑=)()3(204x l x xk k n k k( )。
5、设1326)(247+++=x x x x f 和节点,,2,1,0,2/ ==k k x k 则=],,,[10n x x x f 和=∆07f。
6、5个节点的牛顿-柯特斯求积公式的代数精度为 ，5个节点的求积公式最高代数精度为 。
7、{}∞=0)(k kx ϕ是区间]1,0[上权函数x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x ϕ，则⎰=14)(dx x x ϕ 。
8、给定方程组⎩⎨⎧=+-=-221121b x ax b ax x ，a 为实数，当a 满足 ，且20<<ω时，SOR 迭代法收敛。
9、解初值问题00(,)()y f x y y x y '=⎧⎨=⎩的改进欧拉法⎪⎩⎪⎨⎧++=+=++++)],(),([2),(]0[111]0[1n n n n n n n n n n y x f y x f h y y y x hf y y 是阶方法。
10、设⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=11001a a a a A ，当∈a （ ）时，必有分解式T LL A =，其中L 为下三角阵，当其对角线元素)3,2,1(=i l ii 满足（ ）条件时，这种分解是唯一的。
天大2020年春季考试《数值计算方法》在线作业一.doc
1.解方程组：x1+x2+x3=6;x1+3x2-2x3=1;2x1-2x2+x3=1.x1,x2,x3分别为（）A.1,1,1B.1,1,2C.1,2,3D.3,2,1【参考答案】: C2.若A是n*n阶非奇异阵，则必存在单位下三角阵L和上三角阵U，使A=LU唯一成立A.正确B.错误【参考答案】: B3.设Ax=b,准确解为X*，某一近似解为X，用（）来判断误差A.||AX-b||B.||X-X*||C.bD.||b-AX||【参考答案】: B4.已知多项式P(x)，过点（0,0）（2,8）（4,64）（11,1331）（15,3375），它的三阶差商为常数1，一阶二阶差商均不是0，那么P(x)是（）A.二次多项式B.不超过二次的多项式C.三次多项式D.四次多项式【参考答案】: C5.三次样条函数是在各个子区间上的（）次多项式A.2B.3C.4D.5【参考答案】: B6.若误差限为0.5×10^(-5)，那么近似数0.003400有5位有效数字A.正确B.错误【参考答案】: B7.A.AB.BC.CD.D【参考答案】: C8.区间[a,b]上的三次样条函数是一个次数不超过三次的多项式A.正确B.错误【参考答案】: B9.数值求积公式中的Simpson公式的代数精度为A.0B.1C.2D.3【参考答案】: D10.如果插值结点相同，在满足相同插值条件下所有的插值多项式是等价的。
A.正确B.错误【参考答案】: A11.A.AB.BC.CD.D【参考答案】: A12.下列说法错误的是（）A.如果一个近似数的每一位都是有效数字，则称该近似数为有效数B.凡是经“四舍五入”得到的近似数都是有效数C.数值方法的稳定性是指初始数据的扰动对计算结果的影响D.病态问题是由数学问题本身的性质决定的，与数值方法有关【参考答案】: B13.A.AB.BC.CD.D【参考答案】: D14.双点割线法是用过两点的割线与x轴交点的横坐标作为X0的新近似值A.正确B.错误【参考答案】: A15.设求方程f(x)=0的根的牛顿法收敛，则它具有（）收敛A.超线性B.平方C.线性D.三次【参考答案】: C16.经过A(0,1),B(1,2),C(2,3)的插值多项式P(x)=()A.xB.x1C.2x1D.x^21【参考答案】: D17.A.AB.BC.CD.D【参考答案】: A18.当x=1,-1,2时，f(x)=0,-3,4,则f(x)的二次插值多项式为5x^2/6+2x/3-7/3A.正确B.错误【参考答案】: A19.拉格朗日插值基函数在节点上的取值是（0或1?）A.正确B.错误【参考答案】: A20.x1=1,f(x1)=4;x2=2,f(x2)=1;x3=4,f(x3)=6.则二阶差商f(x1,x2,x3)=()A.5/2B.-3C.11/6D.2/5【参考答案】: C21.含有n+1个节点的插值型数值积分公式的代数精度至少为（）A.1B.nC.n1D.2【参考答案】: B22.若线性方程组的系数矩阵A的条件数Cond（A）=1，则称此方程为良态方程组A.正确B.错误【参考答案】: B23.数值计算中主要研究的误差有（）A.模型误差和观测误差B.截断误差和方法误差C.相对误差和绝对误差 D.模型误差和方法误差【参考答案】: C24.下列求积公式中用到外推技术的是（）A.梯形公式B.复合抛物线公式C.龙贝格公式D.高斯型求积公式【参考答案】: B25.下列说法不正确的是（）A.二分法不能用于求函数f(x)=0的复根B.方程求根的迭代解法的迭代函数为?f(x)，则迭代收敛的充分条件是?f(x)1C.用高斯消元法求解线性方程组AX＝B时，在没有舍入误差的情况下得到的都是精确解D.如果插值节点相同，在满足插值条件下用不同方法建立的插值公式是等价的【参考答案】: B26.在区间[xi-1,xi]上作线性插值，在图形上即为把两点用线段相连，n条线段组成折线，该折线对应的函数称为（）A.牛顿插值函数B.分段线性插值函数C.三次样条插值函数D.拉格朗日插值函数【参考答案】: B27.A.AB.BC.CD.D【参考答案】: C28.牛顿插值多项式的优点是在计算时，高一级的插值多项式可利用前一次插值的结果A.正确B.错误【参考答案】: A29.f(x)=x^2+1,则f[1,2,3,4]=()A.4B.3C.1D.0【参考答案】: D30.A.AB.BC.CD.D 【参考答案】: C31.A.AB.BC.CD.D 【参考答案】: D32.有数表x=0,0.5,1,1.5,2,2.5;f(x)=-2,-1.75,-1,0.25,2,4.25.所确定的插值多项式次数是（）A.二次B.三次C.四次D.五次【参考答案】: A33.若方阵A的谱半径小于1，则解方程组Ax=b的Jacobi迭代法收敛A.正确B.错误【参考答案】: B34.A.AB.BC.CD.D【参考答案】: B35.用二分法求方程x3+x-4=0在区间[1,2]内的根，要求精确到第三位小数，需对分（）次A.8B.11C.9D.10【参考答案】: D36.若线性方程组Ax=b的系数矩阵A为严格对角占优矩阵，则解方程组的Jacobi迭代法和Gauss-Seidel迭代法（）A.都收敛B.都发散C.Jacobi迭代法收敛，Gauss-Seidel迭代法发散 D.Jacobi迭代法发散，Gauss-Seidel迭代法收敛【参考答案】: A37.用简单迭代法求方程f(x)=0的实根，把方程f(x)=0 表示成x=j(x),则f(x)的根是（）A.y=j(x)与x轴交点的横坐标B.y=j(x)与y=x交点的横坐标C.y=x与x 轴的交点的横坐标D.y=j(x)与y=x交点【参考答案】: B38.A.AB.BC.CD.D【参考答案】: C39.A.AB.BC.CD.D【参考答案】: B40.一阶均差f(x0,x1)=(f(x0)-f(x1))/(x0-x1)A.正确B.错误【参考答案】: A。
数值计算大作业
数值计算大作业题目一、非线性方程求根1.题目假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。
（1）如果令()N t 表示在t 时刻的人口数目，β表示固定的人口出生率，则人口数目满足微分方程()()dN t N t dt β=，此方程的解为0()=tN t N e β； （2）如果允许移民移入且速率为恒定的v ，则微分方程变成()()dN t N t vdt β=+， 此方程的解为0()=+(1)t t vN t N e e βββ-；假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率β，精确到410-；且通过这个数值来预测第二年年末的人口数，假设移民速度v 保持不变。
4350001564000=1000000(1)e e βββ+-2.数学原理采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(=x f ，如果)(x f 是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。
设已知方程0)(=x f 有近似根k x (假定0)(≠'x f )，将函数)(x f 在点k x进行泰勒展开，有.))(()()(⋅⋅⋅+-'+≈k k k x x x f x f x f于是方程0)(=x f 可近似地表示为))(()(=-'+k k x x x f x f这是个线性方程，记其根为1k x +，则1k x +的计算公式为)()(1k k k k x f x f x x '-==+，,,2,1,0⋅⋅⋅=k这就是牛顿迭代法，简称牛顿法。
3.程序设计作出函数的图像，大概估计出根的位置fplot('1000*exp(x)+(435*x)*(exp(x)-1)-1564',[0 3]);grid大概估计出初始值x=0.5function [p1,err,k,y]=newton(f,df,p0,delta,max1) % f 是非线性系数 % df 是f 的微商 % p0是初始值% dalta 是给定允许误差 % max1是迭代的最大次数 % p1是牛顿法求得的方程近似解 % err 是p0误差估计 % k 是迭代次数 p0,feval('f',p0) for k=1:max1p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1;p1,err,k,y=feval('f',p1) if(err<delta)|(y==0), break,endp1,err,k,y=feval('f',p1) endfunction y=f(x)y=1000000*exp(x)+435000*(exp(x)-1)/x-1564000; function y=df(x)y=1000000*exp(x)+435000*(exp(x)/x-(exp(x)-1)/x^2);4.结果分析与讨论newton('f','df',1.2,10^(-4),10) 运行后得出结果 p0 =0.5000p1 =0.1679 err =0.3321 k =1 y =9.2415e+004 p1 =0.1031 err =0.0648 k =2 y =2.7701e+003 p1 =0.1010 err =0.0021 k =3 y =2.6953p1 =0.1010 err =2.0129e-006 k =4 y = 2.5576e-006 ans =0.1010运算后的结果为1010.0=β，通过这个数值来预测第二年年末的人口数，0.10100.1010435000f(t)=1000000(1)0.1010t te e +-t=2时候对于f ()2187945.865x =实践表明，当初始值难以确定时，迭代法就不一定收敛了，因此要根据问题实际背景或者二分法先得一个较好的初始值，然后再进行迭代；再者迭代函数选择不合适的话，采用不动点迭代法也有可能出现不收敛的情况；因此我采用的是牛顿法。
数值计算方法习题.doc
第一章 绪论1．把下列各数按四舍五入规则舍入为有3位小数的近似数，并写出近似数的绝对误差和相对误差，指出近似数有几位有效数字： 93.1822 4.32250 15.9477 17.3675 2．按四舍五入原则，将下列数舍成五位有效数字：816.9567 6.000015 17.32250 1.235651 93.18213 0.015236233．设 **,671.3,6716.3x x x 则==有几位有效数字？ 4．若1/4用0.25来表示,问有多少位有效数字? 5．若 1.1062,0.947a b ==是经过舍入后得到的近似值,问:,a b a b +⨯各有几位有效数字?6．设120.9863,0.0062y y ==是经过舍入后作为12,x x 的近似值, 求11y 和21y 的计算值与真值的相对误差限及12y y 和得到真值的相对误差限. 7．设0,x x >的相对误差为δ,求ln x 的绝对误差.8．正方形的边长约为100cm ,应该怎样测量,才能使其面积的误差不超过12cm . 9．设x 的相对误差为a %,求x n 的相对误差.10．计算球的体积,为了使相对误差限为1%,问度量半径R 时允许的相对误差限如何?11．5631.2*=x 是经四舍五入得到的近似值，则其相对误差≤*r e __________ 12． 设 0000073.0,1416.3,1415926.3**=-==x x x x 则称_________误差13．设⎰+=1061dx xx I nn ，设计一个计算10I 的算法，并说明你的算法的合理性。
14．设028Y =，按递推公式1n n Y Y -= （1,2,n = ）， 计算到100Y27.982≈（5位有效数字），试问计算100Y 将有多大误差． 15．求方程25610x x -+=的两个根，使它至少具有4位有效数字27.982≈）．16．当N 充分大时，怎样求121d 1N N x x ++⎰？17．序列{}n y 满足递推关系101n n y y =- （1,2,n = ），若0 1.41y =≈（三位有效数字），计算到10y 时误差有多大？这个计算过程稳定吗？18．计算61)f =1.41≈，利用下列算式计算，哪一个得到的结果最好？，3(3-，99- 19．()ln(f x x =，求(30)f 的值，若开平方用6位函数表，问求对数时误差有多大？若改用另一等价公式ln(ln(x x =- 计算，求对数时误差有多大？第二章 解线性方程组的直接方法1．用高斯消去法解方程组123234011921261x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦ 2．用LU 分解，将上题系数矩阵分解为L 和U 的乘积，L 是对角线元素为1的下三角矩阵，U 是上三角矩阵。
《数值计算方法》试题与答案
习题一1．设x >0相对误差为2%4x 的相对误差。
解：由自变量的误差对函数值引起误差的公式：(())(())'()()()()f x xf x f x x f x f x δδ∆=≈得（1）()f x =11()()*2%1%22x x δδδ≈===；（2）4()f x x =时444()()'()4()4*2%8%x x x x x xδδδ≈===2．设下面各数都是经过四舍五入得到的近似数，即误差不超过最后一位的半个单位，试指出他们各有几位有效数字。
（1）12.1x =；（2）12.10x =；（3）12.100x =。
解：由教材9P 关于1212.m nx a a a bb b =±型数的有效数字的结论，易得上面三个数的有效数字位数分别为：3，4，53．用十进制四位浮点数计算 （1）31.97+2.456+0.1352； （2）31.97+（2.456+0.1352）哪个较精确？解：（1）31.97+2.456+0.1352 ≈21((0.3197100.245610)0.1352)fl fl ⨯+⨯+ =2(0.3443100.1352)fl ⨯+=0.3457210⨯（2）31.97+（2.456+0.1352）21(0.319710(0.245610))fl fl ≈⨯+⨯ = 21(0.3197100.259110)fl ⨯+⨯ =0.3456210⨯易见31.97+2.456+0.1352=0.345612210⨯，故（2）的计算结果较精确。
4．计算正方形面积时，若要求面积的允许相对误差为1%，测量边长所允许的相对误差限为多少？ 解：设该正方形的边长为x ，面积为2()f x x =，由(())(())'()()()()f x xf x f x x f x f x δδ∆=≈解得(())()()'()f x f x x xf x δδ≈=2(())(())22f x x f x x xδδ==0.5%5．下面计算y 的公式哪个算得准确些？为什么？（1）已知1x <<，（A ）11121xy x x-=-++，（B ）22(12)(1)x y x x =++； （2）已知1x>>，（A ）y=，（B ）y = （3）已知1x <<，（A ）22sin x y x =，（B ）1cos2xy x-=；（4）（A）9y =-（B ）y =解：当两个同（异）号相近数相减（加）时，相对误差可能很大，会严重丧失有效数字；当两个数相乘（除）时，大因子（小除数）可能使积（商）的绝对值误差增大许多。
数值计算大作业
课程设计课程名称：设计题目：学号：姓名：完成时间：题目一：非线性方程求根 一 摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。
本实验通过使用常用的求解方法二分法和Newton 法及改进的Newton 法处理几个题目，分析并总结不同方法处理问题的优缺点。
观察迭代次数，收敛速度及初值选取对迭代的影响。
用Newton 法计算下列方程（1） 310x x --= ， 初值分别为01x =，00.45x =，00.65x =； （2） 32943892940x x x +-+= 其三个根分别为1,3,98-。
当选择初值02x =时给出结果并分析现象,当6510ε-=⨯，迭代停止。
解：1)采用MATLAB 进行计算；首先定义了Newton 法：function kk=newton(f,df,x0,tol,N)% Newton Method （牛顿法）% The first parameter f is a external function with respect to viable x.（第一个参数也就是本题所用的函数f ）% The second parameter df is the first order diffential function of fx.（第二个参数也就是本体所用函数f 的导数方程df ） % x0 is initial iteration point(初值). % tol is the tolerance of the loop （精度）.% N is the maximum number of iterations （循环上限）. x=x0;f0=eval(f);df0=eval(df); n=0;disp(' [ n xn xn+1 fn+1 ]'); while n<=N x1=x0-f0/df0; x=x1; f1=eval(f); X=[n,x0,x1,f1]; disp(X);if abs(x0-x1)<tolfprintf('The procedure was successful.') kk=X; return else n=n+1; x0=x1;f0=f1;endendif n==N+1fprintf('the method failed after N iterations. '),kk=0;End我们把Newton法存为.m格式的文件；之后我们运行程序：clear;clc;syms xf=x^3-x-1;df=diff(f,x);x=newton(f,df,1,0.0001,50);x会得到一下结果[ n xn xn+1 fn+1 ]0 1.0000 1.5000 0.87501.0000 1.5000 1.0625 -0.86302.0000 1.0625 1.4940 0.8408到第50次迭代时候会出现该问题：47.0000 1.4898 1.0814 -0.816748.0000 1.0814 1.4898 0.816749.0000 1.4898 1.0814 -0.816750.0000 1.0814 1.4898 0.8167the method failed after N iterations.x =0；同样测试x0=0.45、0.65得不出结果，判断出初值离真值太远，所以我们采用牛顿下山法进行计算迭代：我们定义了其中的f函数和df函数，并且分别存为.m格式的文件，其代码如下：f:function y=f(x)y=x^3-x-1;df:function y=df(x)y=3*x^2-1;之后我们定义newton下山法同时也存为.m的程序：function [x,i]=downnewton(f,df,x0,tol)k=0;i=1;disp(' [ n xn xn+1 fn+1 ]'); while(k==0)fx=feval('f',x0);dfx=feval('df',x0);t=0;u=1;while(t==0)dx=-fx/dfx;x1=x0+u*dx;fx1=feval('f',x1);fx0=feval('f',x0);if(abs(fx1)>abs(fx0));u=u/2;elset=1;endendX=[i,x0,x1,fx1];disp(X);if(abs(fx1)<tol)k=1;elsex0=x1;i=i+1;endendx=x1;i=i;end之后带入x0=0.45;downnewton('f','df',0.45,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.4500 -0.4155 -0.65622.0000 -0.4155 -0.5857 -0.61523.0000 -0.5857 -0.5754 -0.61514.0000 -0.5754 -0.5782 -0.61515.0000 -0.5782 -0.5773 -0.61516.0000 -0.5773 -0.5774 -0.61517.0000 -0.5774 -0.5773 -0.61518.0000 -0.5773 -0.5774 -0.61519.0000 -0.5774 -0.5774 -0.615110.0000 -0.5774 -0.5774 -0.615111.0000 -0.5774 1.3131 -0.049012.0000 1.3131 1.3248 0.000513.0000 1.3248 1.3247 0.0000ans =1.3247带入x0=0.6;downnewton('f','df',0.6,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.6000 1.1406 -0.65662.0000 1.1406 1.3668 0.18663.0000 1.3668 1.3263 0.00674.0000 1.3263 1.3247 0.00005.0000 1.3247 1.3247 0.0000ans =1.3247带入x0=1;downnewton('f','df',1,10^(-6))[ n xn xn+1 fn+1 ]1.0000 1.0000 1.5000 0.87502.0000 1.5000 1.3478 0.10073.0000 1.3478 1.3252 0.00214.0000 1.3252 1.3247 0.0000ans =1.32472)同样采用Newton下山法：重新定义f、df:f:function y=f(x)y=x^3+94*x^2-389*x+294;df:function y=df(x)y=3*x^2+188*x-389;再带入初值x0=2;downnewton('f','df',2,5*10^(-6))[ n xn xn+1 fn+1 ]1 2 -98 0ans =-98得出x=-98；分析：先画出该函数的图像；x=(-100:.1:100);ezplot('x^3+94*x^2-389*x+294',[-100 100]) 得出该图像如图：-100-80-60-40-20020406080024681012141618x 105xx 3+94 x 2-389 x+294根据牛顿法的几何解释，在x0=2的点做切线，与y 相交，交点的横坐标值为x=-98则结束了该现象。
北理工数值分析大作业
数值分析上机作业第 1 章1.1计算积分，n=9。
（要求计算结果具有6位有效数字）程序：n=1:19;I=zeros(1,19);I(19)=1/2*((exp(-1)/20)+(1/20));I(18)=1/2*((exp(-1)/19)+(1/19));for i=2:10I(19-i)=1/(20-i)*(1-I(20-i));endformat longdisp(I(1:19))结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：当计算到数列的第10项时，所得的结果即为n=9时的准确积分值。
取6位有效数字可得.1.2分别将区间[-10.10]分为100,200,400等份，利用mesh或surf命令画出二元函数z=的三维图形。
程序：>> x = -10:0.1:10;y = -10:0.1:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.1')>> x = -10:0.2:10;y = -10:0.2:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长 0.2')>>x = -10:0.05:10;y = -10:0.05:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.05')结果截图及分析：由图可知，步长越小时，绘得的图形越精确。
数值计算方法计算习题
1.已知ln(2.0)=0.6931;ln(2.2)=0.7885,ln(2.3)=0.8329, 试用线性插值和抛物插值计算.ln2.1的值并估计误差（牛顿插值和拉格朗日插值）2.已知函数y=sinx 的数表如下,分别用前插和后插公式计算sin0.57891的值，并估算误差。
i x0.4 0.5 0.6 0.7 )(i x f0.389420.479430.564640.644223. 已知i x-2 -1 0 1 2 )(i x f42135求)(x f 的二次拟合曲线)(2x p ，并求)0(f '的近似值。
4. 数值积分公式形如⎰'+'++=≈1)1()0()1()0()()(f D f C Bf Af x S dx x xf 试确定参数D C B A ,,,使公式代数精度尽量高；（2）设]1,0[)(4C x f ∈，推导余项公式⎰-=1)()()(x S dx x xf x R ，并估计误差。
5. 已知数值积分公式为：)]()0([)]()0([2)(''20h f f h h f f hdx x f h-++≈⎰λ，试确定积分公式中的参数λ，使其代数精确度尽量高，并指出其代数精确度的次数。
6. 用复化Simpson 公式计算积分()⎰=10sin dx x x I 的近似值，要求误差限为5105.0-⨯。
7. 已知012113,,424x x x ===，给出以这3个点为求积节点在[]0.1上的插值型求积公式。
8. 给出 900,cos ≤≤x x 的函数表，步长)60/1(1='=h ，若函数具有5位有效数字，研究用线性插值求x cos 近似值时的总误差界。
9. 求一个次数不高于4次的多项式)(x P ，使它满足0)0()0(='=P P ，1)1()1(='=P P ，1)2(=P 。
10. 单原子波函数的形式为bxae y -=，试按照最小二乘法决定参数a 和b ，已知数据如下：X 0 1 2 4 y2.010 1.210 0.740 0.45011. 分别用梯形公式和辛普森公式计算下列积分：⎰+1024dx x x。
龙贝格-求积分
数值计算方法大作业注:由朱福利邹素云共同完成龙贝格求积公式【简介】龙贝格求积公式也称为逐次分半加速法。
它是在梯形公式、辛卜生公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。
这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
【算法】对区间[a， b]，令h=b-a构造梯形值序列{T2K}。
T1=h[f(a)+f(b)]/2把区间二等分，每个小区间长度为 h/2=(b-a)/2，于是T2 =T1/2+[h/2]f(a+h/2）把区间四(2)等分，每个小区间长度为h/2 =（b-a）/4，于是T4 =T2/2+[h/2][f（a+h/4）+f（a+3h/4).....................把[a，b] 2等分，分点xi=a+（b-a）/ 2 ·i （i =0，1，2 · · · 2k）每个小区间长度为（b-a）/ 2 .例：I = ∫0(4/1+X) dx解按上述五步计算，此处 f(x)=4/(1+x) a=0 b=1 f(0)=4 f(1)=2由梯形公式得T1=1/2[f(0)+f(1)]=3计算f(1/2)=16/5 用变步长梯形公式得T2=1/2[T1+f(1/2)]=3.1由加速公式得S1=1/3(4T2-T1)=3.133333333求出f(1/4) f(3/4) 进而求得T4=1/2{T2+1/2[f(1/4)+f(3/4)]}=3.131176471S2=1/3(4T4-T2)=3.141568628C1=1/15(16S2-S1)=3.142117648计算f(1/8) f(3/8) f(5/8) f(7/8)进而求得T8=1/2{T4+1/4[f(1/8)+f(3/8)+f(5/8)+f(7/8)]}=3.138988495S4=1/3(4T3-T4)=3.141592503C2=1/15(16S4-S2)=3.141594095R1=1/63(64C2-C1)=3.141585784把区间再二分，重复上述步骤算得T16=3.140941613 S8=3.141592652C4=3.141592662 R2=3.141592640程序步骤1.输入积分上线2.输入积分下线3.输入区间等分数4.输入要求的函数5.计算出所求函数的积分，分别是●复化梯形求积结果●辛普森求积结果●科特斯求积结果●龙贝格求积结果流程图源程序// MyT ask.cpp : Defines the entry point for the console application.//#include "stdafx.h"#include "math.h"#include<stdlib.h>#define IS_INT(x) ( (x) == int( (x) * 10 / 10 ))void fuHuaTiXing(int k, int a, int b);//复化梯形求积里面T1跟T2double f(double x);//原函数sinx / xdouble fuHuaTiXing(int k, double y,int a, int b);//复化梯形求积里面T1跟T2double leiJia(int a, int b, int i);//复化梯形求积里面的累加void simpson(int k);//k代表要计算的simpson的个数void cotes(int k);//k >=3void romberg(int k );double TTT[512] = {0.0};double TT[128] = {0.0};double S[64] = {0.0};double C[32] = {0.0};double R[16] = {0.0};int flag;int main(int argc, char* argv[]){int a,b,k;//a 积分下线b 积分上限k 等分数//k最大支持256等分int _k;double result;int loop_k;int g_quit = 0;int g_quit1 = 0;//欢迎界面printf("\t\t**************************************** ************\n");printf("\t\t**************************************** ************\n");printf(" 欢迎使用RomBerg方法计算积分\n");printf(" 本程序支持最小4等分最大256等分\n");printf("\t\t****************************************************\n");printf("\t\t**************************************** ************\n");printf("想继续运行吗，继续请按1，退出请按0\n");scanf("%d",&g_quit);if(1 == g_quit){while(1){printf("please input 3 numbers whitch means three constent.\n");printf("please input the first number whitch means 积分下限\n");getchar();scanf("%d",&a);printf("please input the second number whitch means 积分上限\n");getchar();scanf("%d",&b);printf("please input the third number whitch means 等分数\n\t\t\t\t\t注意:支持4等分到256等分\n \t\t\t\t\t注意:一定是2的整数幂\n");scanf("%d",&k); //if(4 == (int)k){_k = 2;}else if(8 == (int)k) {_k = 3;}else if(16 == (int)k) {_k = 4;}else if(32 == (int)k) {_k = 5;}else if(64 == (int)k) {_k = 6;}else if(128 == (int)k){_k = 7;}else if(256 == (int)k){_k = 8;}else{printf("您输入的等分区间不是2的整数幂,请重新输入!\n\n");getchar();getchar();getchar();continue;}aaaaa:printf("please select the function of f(x). \n");printf("if you want to chose function 1( pow(x,2)*sin(x) ) please input 1;\n");printf("if you want to chose function 2 ( sin(x)*cos(x) ) please input 2;\n");printf("if you want to chose function 3 ( sin(x)*sin(x)*x ) please input 3;\n");printf("if you want to chose function 4 ( sin(x)/x ) please input 4.\n");scanf("%d",&flag);if(flag==1 || flag==2 || flag==3 || flag==4 ){//null}else{printf("You have inputed wrong num! \nPlease try again!\n\n\n");goto aaaaa;}fuHuaTiXing(_k,a,b);simpson(_k);cotes(_k);romberg(_k);printf("output fuHuaTiXing(复化梯形) result.\n");for(loop_k=0;loop_k<=_k;loop_k++){printf("TT%d=%.14f\n",(int)pow(2,loop_k),TT[loop_k]);}printf("output simpson result.\n");for(loop_k=0;loop_k<=_k-1;loop_k++){printf("S%d= %.14f\n",(int)pow(2,loop_k),S[loop_k]);}printf("output cotes result.\n");for(loop_k=0;loop_k<=_k-2;loop_k++){printf("C%d = %.14f \n",(int)pow(2,loop_k),S[loop_k]);}printf("output romberg result.\n");for(loop_k=0;loop_k<=_k-3;loop_k++){printf("R%d= %.14f\n",(int)pow(2,loop_k),R[loop_k]);}printf("\n\n\n");printf("继续计算请按1，退出请按0\n");scanf("%d",&g_quit1);if(0 == g_quit1){exit(0);}}}else{exit(0);}return 0;}double f(double x)//原函数sinx/x {double result_f;switch(flag){case 1:result_f = pow(x,2)*sin(x);return (result_f);break;case 2:result_f = sin(x)*cos(x);return (result_f);break;case 3:result_f = sin(x)*sin(x)*x;return (result_f);break;case 4:result_f = sin(x);if(x != 0){result_f = result_f / x;return (result_f);}else return 1.0;break;default :break;}}void fuHuaTiXing(int k, int a, int b)//复化梯形求积里面T1跟T2{int loop_i;double result_b_a;result_b_a = double(b - a) / 2;for(loop_i = 0; loop_i<=k; loop_i++)//k >=3{int temp;if(loop_i == 0)//k means zhishu{TT[0] = result_b_a * ( f(a) + f(b) );//代表T1}else{temp = pow(2,loop_i);TT[loop_i] = TT[loop_i - 1] / 2 + ( (double)(b - a) / (temp ) ) * leiJia(a,b,loop_i);}}for(loop_i=0;loop_i<10;loop_i++){TTT[loop_i] = TT[loop_i];//把局部数组转存到全局数组中方便下面的访问}}double leiJia(int a, int b, int i)//复化梯形求积里面的累加{int loop_m ;int loop_i;double result_leijia = 0.0;int temp;double temp1,temp2;switch(i){case 1:loop_m = (int)pow(2,(i-1));break;case 2:loop_m = (int)pow(2,(i-1));break;case 3:loop_m = (int)pow(2,(i-1));break;case 4:loop_m = (int)pow(2,(i-1));break;case 5:loop_m = (int)pow(2,(i-1));break;case 6:loop_m = (int)pow(2,(i-1));break;default:break;}for(loop_i = 1; loop_i<=loop_m; loop_i++){temp = pow(2,i);temp1 = (double)(a + (2 * loop_i - 1) * ( (double)(b - a) / temp ) );//temp2 = f( temp1 );result_leijia = result_leijia + f( temp1 ); //第一次temp1=0.5}return (result_leijia);}//求simpsonvoid simpson(int k)//k代表要计算的ximpson的个数{int loop_i = 0, loop_j = 0;for(loop_i =0;loop_i<=k-1;loop_i++)S[loop_i] = ( 4 * TT[loop_i + 1] -TT[loop_i] ) / 3;}void cotes(int k){int loop_i;for(loop_i=0;loop_i<=k-2;loop_i++)C[loop_i] = ( 16 * S[loop_i+1] - S[loop_i] ) / 15;}void romberg(int k ){int loop_i;for(loop_i=0;loop_i<=k-3;loop_i++)R[loop_i] = ( pow(4,3) * C[loop_i +1] - C[loop_i] ) / (pow(4,3) - 1) ;}。
数值计算方法丁丽娟课后习题答案
数值计算方法丁丽娟课后习题答案【篇一：北京理工大学数值计算方法大作业数值实验1】）书p14/4分别将区间[?10,10]分为100,200,400等份，利用mesh或surf命令画出二元函数的三维图形。
z=|??|+ ??+?? +??++??【matlab求解】[x,y]=meshgrid(-10:0.1:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.05:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.025:10); a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);（二）书p7/1.3.2数值计算的稳定性（i）取= ??c语言程序—不稳定解 +=ln1.2，按公式=?? （n=1,2，…） #includestdio.h#includeconio.h#includemath.hvoid main(){float m=log(6.0)-log(5.0),n;int i;i=1;printf(y[0]=%-20f,m); while(i20){n=1/i-5*m;printf(y[%d]=%-20f,i,n);m=n;i++;if (i%3==0) printf(\n); }getch();}（ii） c语言程序—稳定解≈??[ ??+?? +?? ??+??按公式 =??(??)#includestdio.h#includeconio.h#includemath.hvoid main(){float m=(1/105.0+1/126.0)/2,n; k=n,n-1,n-2，…）（【篇二：北京理工大学数值计算方法大作业数值实验4】 p260/1考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为= ?????????????????? ?? ??????????= ?????????????? ??曲线关于原点对称，取a=1，参数s的变化范围[-5,5]，容许误差限分别是，,和。
数值计算方法第1章作业
第一章作业第一题问题叙述：构造算法并编程序精确计算二次方程的根。
●设a≠0，b2-4ac>0，且有方程ax2+bx+c=0●包括b2≈b2-4ac的情况（a=c=1，b=±1000000.000001）问题分析：对于普通的二次求根公式：x1,2=−b±√b2−4ac2a当b2>>4ac时，分子可能非常小。
由于计算机中的算术运算存在减性抵消的现象，即两个几乎相等的浮点数相减时会引起舍入误差，所以在这种极端条件下用这个公式就会带来很大的误差。
解决方法：1.使用双精度2.使用变换公式x1,2=−2cb±√b2−4ac3.先利用原公式计算较大的根（即分子不会引起减性抵消），再利用公式：x1x2=c a计算较小的根。
问题解决：1.使用双精度：%This program uses double precision to solve the equation of two degree%And make a comparision to the single precisionclear;clc;[a,b,c]=textread('data.txt','%n%n%n'); %read the numbers from data.txtdelta=b*b-4*a*c;x1=(-b+sqrt(delta))/(2*a);x2=(-b-sqrt(delta))/(2*a); %double precisiona=single(a);b=single(b);c=single(c);%use single precisiondelta=single(b*b-4*a*c);x11=single((-b+sqrt(delta))/(2*a));x12=single((-b-sqrt(delta))/(2*a));fid=fopen('out.txt','w');fprintf(fid,'%g %g %g %g',[x1 x2 x11 x12]);fclose(fid);%write the result into the out.txt下面是计算结果：结论：1.在一般情况下，即没有出现b2≈b2−4ac时，无论是单精度还是双精度下均可以得出正确答案。
数值计算方法大作业
目录第一章非线性方程求根 (3)1.1迭代法 (3)1.2牛顿法 (4)1.3弦截法 (5)1.4二分法 (6)第二章插值 (7)2.1线性插值 (7)2.2二次插值 (8)2.3拉格朗日插值 (9)2.4分段线性插值 (10)2.5分段二次插值 (11)第三章数值积分 (13)3.1复化矩形积分法 (13)3.2复化梯形积分法 (14)3.3辛普森积分法 (15)3.4变步长梯形积分法 (16)第四章线性方程组数值法 (17)4.1约当消去法 (17)4.2高斯消去法 (18)4.3三角分解法 (20)4.4雅可比迭代法 (21)4.5高斯—赛德尔迭代法 (23)第五章常积分方程数值法 (25)5.1显示欧拉公式法 (25)5.2欧拉公式预测校正法 (26)5.3改进欧拉公式法 (27)5.4四阶龙格—库塔法 (28)数值计算方法第一章非线性方程求根1.1迭代法程序代码：Private Sub Command1_Click()x0 = Val(InputBox("请输入初始值x0"))ep = Val(InputBox(请输入误差限ep))f = 0While f = 0X1 = (Exp(2 * x0) - x0) / 5If Abs(X1 - x0) < ep ThenPrint X1f = 1Elsex0 = X1End IfWendEnd Sub例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10）1.2牛顿法程序代码：Private Sub Command1_Click()b = Val(InputBox("请输入被开方数x0"))ep = Val(InputBox(请输入误差限ep))f = 0While f = 0X1 = x0 - (x0 ^ 2 - b) / (2 * b)If Abs(X1 - x0) < ep ThenPrint X1f = 1Elsex0 = X1End IfWendEnd Sub例：求56的值。
计算方法
SHANGHAI JIAO TONG UNIVERSITY论文题目：计算方法大作业多方法求解数值积分****: ***学生学号: ************专业: 机械工程****: ***学院(系): 机械与动力工程学院多方法求解数值积分具体题目要求：用不同数值方法计算积分49xdx=-⎰(1) 取不同的步长h，分别用复合梯形及复合辛普森公式计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？(2) 用龙贝格求积计算完成问题(1);(3) 用自适应辛普森积分，使其精度达到410-。
一．设计目的由积分学基本理论，定积分可由Newton Leibniz-公式计算，但是对于一些无法找到原函数的函数（如2xe-等）不能通过牛顿—莱布尼兹公式计算，就必须得另寻它法。
因此需要我们能够熟练地应用常用的数值积分计算方法（如机械求积、Newton Cotes-公式等）并掌握结合数值计算软件（Matlab、Lingo 等）及计算机高级语言()c java、进行对应算法实现的技能。
熟练数学软件求解数学问题，掌握各种数学问题的求解方法。
本设计主要是通过多种复合求积公式求解积分，主要包括复化梯度法、复化辛普森法、龙贝格以及自适应辛普森法等求解方法，利用C#语言编写相对应的算法进行求解，并用Matlab作图分析，大大地提高了解题的速度。
二．积分算法1．复合梯形公式、复合辛普森公式算法当积分区间[a,b]的长度较大，而节点数n+1固定时，直接使用牛顿-柯特斯公式的余项将会较大，而如果增加节点个数，即n+1增加时，公式的舍入误差又很难得到控制。
为了提高公式的精度，又使算法简单易行，往往使用复合方法，即将积分区间[a,b]分成若干个子区间，然后在每个小区间上的积分的近似值相加将定积分∫f(x)dxba的积分区间[a,b]分割n等分，各节点为x k=a+kℎ,k=0,1⋯n ℎ=b−an在子区间[x k ,x k+1] (k =0,1⋯n −1 )上使用牛顿-柯特斯公式将[x k ,x k+1]分割为l 等份，步长为ℎl ，节点为x k ,x k +ℎl ,x k +2ℎl ,⋯,x k +lℎl=x k+1记为x k ,x k+1l,x k+2l,⋯,xk+ll=x k+1在[x k ,x k+1]上作f(x)的l 阶牛顿-柯特斯公式∫f(x)dx x k+1x k≈I i (k )=(x k+1−x k )∑C i (l )f (xk+i l)li=0=ℎ∑C i (l )f(xk+i l)li=0由积分的区间可加性，可得∫f(x)dx b a=∑∫f(x)dx x k+1x kn−1k=0≈∑I i (k )n−1k=0=ℎ∑∑C i (l )f(xk+i l)li=0=I n n−1k=0l =1时，可得复合梯形求积公式∫f(x)dx ba≈T n =ℎ∑12n−1k=0[f (x k )+f (x k+1)]即复合梯形公式T n =b −a2n[f (a )+2∑f (x k )n−1k=1+f (b )] 2()12b a E h f η-''=-[,]a b η∈ l =2时，可得复合辛普森求积公式∫f(x)dx ba≈S n =ℎ∑16n−1k=0[f (x k )+4f (x k+12)+f (x k+1)]即复合辛普森公式S n =b −a6n [f (a )+4∑f (x k+12)n−1k=0+2∑f (x k )n−1k=1+f (b )]4(4)()1802b a h E f η-⎛⎫=- ⎪⎝⎭[,]a b η∈算法流程图如下：2．龙贝格算法这里，梯形公式显得算法简单，具有如下递推关系121021()22n n n i i h T T f x -+==+∑因此，很容易实现从低阶的计算结果推算出高阶的近似值，而只需要花费较少的附加函数计算。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
题目利用数值计算方法求取基尼系数姓名与学号指导教师年级与专业所在学院一、问题综述：基尼系数（Gini coefficient），是20世纪初意大利学者科拉多·吉尼根据劳伦茨曲线所定义的判断收入分配公平程度的指标。
是比例数值，在0和1之间。
基尼指数（Gini index）是指基尼系数乘100倍作百分比表示。
在民众收入中，如基尼系数最大为“1”，最小等于“0”。
前者表示居民之间的收入分配绝对不平均（即所有收入都集中在一个人手里，其余的国民没有收入），而后者则表示居民之间的收入分配绝对平均，即人与人之间收入绝对平等，但这两种情况只出现在理论上；因此，基尼系数的实际数值只能介于0～1之间，基尼系数越小收入分配越平均，基尼系数越大收入分配越不平均。
设右图中的实际收入分配曲线（红线）和收入分配绝对平等线（绿线）之间的面积为A，和收入分配绝对不平等线（蓝线）之间的面积为B，则表示收入与人口之间的比例的基尼系数为AA+B。
如果A为零，即基尼系数为0，表示收入分配完全平等（红线和绿线重叠）；如果B为零，则系数为1，收入分配绝对不平等（红线和蓝线重叠）。
该系数可在0和1之间取任何值。
实际上，一般国家的收入分配，既不是完全平等，也不是完全不平等，而是在两者之间，劳伦茨曲线为一条凸向横轴的曲线。
收入分配越趋向平等，劳伦茨曲线的弧度越小（斜度越倾向45度），基尼系数也越小；反之，收入分配越趋向不平等，劳伦茨曲线的弧度越大，那么基尼系数也越大。
基尼系数的调节需要国家通过财政政策进行国民收入的二次分配，例如对民众的财政公共服务支出和税收等，从而让收入均等化，令基尼系数缩小。
基尼系数由于给出了反映居民之间贫富差异程度的数量界线，可以较客观、直观地反映和监测居民之间的贫富差距，预报、预警和防止居民之间出现贫富两极分化。
因此得到世界各国的广泛认同和普遍采用。
联合国有关组织规定：●若低于0.2表示收入平均；●0.2-0.3表示相对平均；●0.3-0.4表示相对合理；●0.4-0.5表示收入差距大；●0.6以上表示收入差距悬殊。
2013年1月18日，中国国家统计局一次性公布了自2003年以来十年的全国基尼系数。
大陆统计局局长马建堂称，按照国际新的统计口径，大陆居民收入的基尼系数，2003年是0.479，2004年是0.473，2005年为0.485，2006年为0.487，2007年为0.484，2008年为0.491，2009年为0.490，2010年为0.481，2011年为0.477，到2012年的数据是0.474，为2005年以来最低水平，而自2008年起，基尼系数也在逐年下降。
而此前西南财大调查数据显示，中国的2012年的基尼系数为0.61，但无论是民间统计的数据还是官方统计的数据，结果都遭到学术界质疑，仍具有争议性。
本文将根据网络上国家统计局的数据，利用上面给出的公式来计算我国从2002年以来的城镇居民基尼系数，并将计算出的数据与现有数据进行比较。
全球基尼系数二、问题分析：想要通过拟合曲线再积分求取基尼系数，需要先求取劳伦茨曲线。
劳伦茨曲线是1905年由经济学家马克斯·劳伦茨所提出的表示收入分配的曲线。
在经济学里,劳伦兹曲线是在过往财富分配数据上建立的累积分布函数所对应的曲线,它通过变量y%的值来反映各项分配的比例。
它经常被用来描述收入的分配情况，即以x%代表一部分(收入相似)家庭占整个社会家庭的比例，以y%代表该部分家庭的收入占整个社会收入的比例。
该曲线也可用来描述社会资本的分配情况。
在这些应用当中，经济学家经常把它用来衡量社会(主要指社会收入)是否公平。
所利用的概率密度分布函数为：本文中使用国家统计局网站上面的数据，先采用拟合得方法来得到劳伦茨曲线，然后再用积分的方法来计算我国城镇居民的基尼系数。
并分别使用多种方法进行计算，例如直接用求得的数据进行插值积分来计算基尼系数，比较不同，再和上文中所提到的国家统计局统计出的基尼系数进行比较，进行进一步的分析。
三、问题解决首先在国家统计局的网站上面找到了按收入等级分城镇居民家庭基本情况——城镇居民平均每人全部年收入。
列举数据表如下（困难户是包含在最低收入户中的群体）：表1.2007-2012年城镇居民平均每人全部年收入表2.2002-2006年城镇居民平均每人全部年收入统计局的数据将所有城镇居民按收入划分为八个部分，按照劳伦茨曲线的定义，我们对每年的数据进行分析都可以得到七个数据点。
通过这七个数据点再加上固有的（0,0），（1,1）点，每年共9个点，就是我们需要的数据。
数据处理方法：从低到高算起，拿2012年的数据为例，对于城镇居民困难户（5%），其收入占总收入的比例为：7520.90∗5%=1.39%26959.00∗100%如果令第一个点为固定点（0,0），则由上面的结果可以求得第二个点为：（0.0139,0.05）。
同样，对于城镇居民最低收入户(10%)，其收入占总收入的比例为：9209.50∗10%=3.42%26959.00∗100%由于困难户是包含在最低收入户中的，所以不用进行累加，可直接得到第三个点：（0.0342,0.10）。
之后，对于城镇居民较低收入户(10%)，其收入占总收入的比例为：13724.70∗10%=5.09%26959.00∗100%由于与之前不再具有包含关系，所以需要与之前的计算结果进行累加，可得到第四个点：（0.0851,0.20）。
之后的计算方式以此类推，此处不再赘述。
以此种方法对以上数据处理后，可得到如下结果：表3.2007-2012年数据处理结果表4.2002-2006年数据处理结果1.最小二乘法拟合最小二乘法是拟合近似曲线，使偏差平方和最小。
一般设：φ(x)=a0φ0(x)+a1φ1(x)+⋯+a nφn(x)mS(a0,a1,……,a n)=∑[a0φ0(x)+a1φ1(x)+⋯+a nφn(x)−y i]2i=1求S的最小值，即∂S=0⇒a0(φk,φ0)+a1(φk,φ1)+⋯+a n(φk,φn)=(φk,f)k其中k=0,1…n对于代数多项式拟合，可以得到：源程序zxecf.m：function[a]=zxecf(x,y,n)% input: x is independent variable matrix% y is dependent variable matrix% n is the order of the fitting% output:a is the polynomial coefficient matrix of fittingm=length(x);C=zeros(m,n+1);for i=1:n+1C(:,i)=x'.^(i-1);endA=C'*C;b=C'*y';a=gauss(A,b); %Gaussian eliminationfor i=1:n+1t(i)=a(n-i+2);enda=t;上面代码中用到的高斯消去法的代码如下：源程序gauss.m：function [x]=gauss(a,b)% input the A,B of AX=B% output the X of AX=Bn=length(a);ab = [a b]; % augmented matrixfor k=1:n-1[~,index] = max(abs(ab(k:n,k))); %gain the location of max numberindex=index+k-1; % The biggest column main entry number of rowsif index ~= ktemp = ab(index,:);ab(index,:) = ab(k,:);ab(k,:) = temp;end%exchange the column main entry to current rowfor j=k+1:nfactor=ab(j,k)/ab(k,k);ab(j,k:n+1) = ab(j,k:n+1) - factor.*ab(k,k:n+1);endend%elimination processx(n)=ab(n,n+1)/ab(n,n);for i=n-1:-1:1sum=ab(i,n+1);for j=i+1:nsum=sum-ab(i,j)*x(j);endx(i)=sum/ab(i,i);endend以2012年数据为例，输入：先试进行四阶拟合，运行源程序a=zxecf(x,y,4)，可得：常数项不为0，可见拟合效果不是很好，继续增大拟合阶数至7阶时，可得：即y=−6.0662x7+20.2785x6−26.4926x5+17.7751x4−6.6249x3+ 1.9313x2+0.1988x此时拟合效果比较好。
拟合效果如下图：x=[0 0.05 0.1 0.2 0.4 0.6 0.8 0.9 1]y=[0 0.0139 0.0342 0.0851 0.2214 0.4034 0.6464 0.8077 1]a =0.6832 -0.9742 1.0406 0.2511 -0.0004a =-6.0662 20.2785 -26.4926 17.7752 -6.6249 1.9313 0.1988 -0.0000对所有年份进行七阶拟合，结果如下表所示：表5.最小二乘法拟合结果2.龙贝格积分龙贝格积分法是通过用若干个积分近似值来推算更精确的新的近似值的方法。
外推公式为：R(j,k)=R(j,k−1)+R(j,k−1)−R(j−1,k−1)4k−1其中第一列：T0(0)=b−a2[f(a)+f(b)]T0(l)=1T0(l−1)+b−al∑f[a+(2i−1)b−al]2l−1i=1,l=1,2,3…外推到k=4即可，因为进一步外推，外推值趋向于0，此时意义不大，所以可以终止进一步的外推。
代码romberg.m：function [i,r,gi]=romberg(a,n)% input:a is the poly coefficient matrix% n is the order of integral% output:i--integral value% r is romberg matrix% gi is the Gini coefficientf=@(x1)polyval(a,x1);x=0;y=1;h=y-x; % h is the intal steplengthr=zeros(n+1,n+1); % initialize the romberg matrixr(1,1)=h*(f(x)+f(y))/2;for j=2:n+1; % calculate the first lineh=h/2;r(j,1)=r(j-1,1)/2+h*sum(f(x+h*(1:2:2^(j-1))));endfor k=2:n+1 % calculate the romberg matrixr(k:n+1,k)=r(k:n+1,k-1)+(r(k:n+1,k-1)-r(k-1:n,k-1))/(4^(k-1)-1);endi=vpa(r(n+1,n+1),5);gi=(0.5-i)/0.5;将上面得到的2012年的系数矩阵输入，可以得到结果：i =0.34238r =0.5001 0 0 0 00.3909 0.3545 0 0 00.3568 0.3455 0.3449 0 00.3462 0.3426 0.3424 0.3424 00.3433 0.3424 0.3424 0.3424 0.3424 gi =0.31524714285717436723643913865089可以得到2012年我国城镇居民的基尼系数为0.3152。