数值分析所有代码

1、实验一：拉格朗日插值多项式命名（源程序.cpp及工作区.dsw）：lagrange问题：0.40.550.800.910.410750.578150.888111.026521.17520用四次拉格朗日多项式求的函数近似值/Lagrange.cpp#include #include #define N 4int checkvalid(double x, int n);void printLag (double x, double y, double varx, int n);double Lagrange(double x, double y, double varx, int n);void

2、 main ()double xN+1 = 0.4, 0.55, 0.8, 0.9, 1;double yN+1 = 0.41075, 0.57815, 0.88811, 1.02652, 1.17520;double varx = 0.5;if (checkvalid(x, N) = 1)printf(nn插值结果: P(%f)=%fn, varx, Lagrange(x, y, varx, N);elseprintf(结点必须互异);getch();int checkvalid (double x, int n)int i,j;for (i = 0; i n; i+)for (j = i

3、+ 1; j n+1; j+)if (xi = xj)/若出现两个相同的结点，返回-1return -1;return 1;double Lagrange (double x, double y, double varx, int n)double fenmu;double fenzi;double result = 0;int i,j;printf(Ln(x) =n);for (i = 0; i n+1; i+)fenmu = 1;for (j = 0; j n+1; j+)if (i != j)fenmu = fenmu * (xi - xj);printf(t%f, yi / fenmu

4、);fenzi = 1;for (j = 0; j n+1; j+)if (i != j)printf(*(x-%f), xj);fenzi = fenzi * (varx - xj);if (i != n)printf(+n);result += yi / fenmu * fenzi;return result;运行及结果显示： 实验二：牛顿插值多项式命名（源程序.cpp及工作区.dsw）：newton_cz问题：0.40.550.800.910.410750.578150.888111.026521.17520用牛顿插值多项式求 /newton_cz.cpp#include#include

5、#include#define N 4int checkvaild(double x,int n)int i,j;for(i=0;in+1;i+) for(j=i+1;jn+1;j+)if(xi=xj)return -1;return 1;void chashang(double x,double y,double fN+1)int i,j,t=0;for(i=0;iN+1;i+)fi0=yi;/f00=y0,f10=y1;f20=y2;f30=y3;f40=y4for(j=1;jN+1;j+)/ 阶数jfor(i=0;iN-j+1;i+)/差商个数ifij=(fi+1j-1-fij-1)/(

6、xt+i+1-xi);/一阶为f01f31;二阶为f02f22 /三阶为f03f13;四阶为f04 t+;double compvalue(double tN+1,double x,double varx)int i,j,r=0;double sum=t00,mN+1=sum,1,1,1,1;printf(itXit F(Xi)t 1阶t 2阶tt3阶t 4阶n);printf(-);for(i=0;iN+1;i+)r=i;printf(x%dt%f ,i,xi);for(j=0;j=i;j+)printf(%f ,trj);r-;printf(n); printf(-);/*/printf(

7、N(x) =n);printf( %fn,t00);for(i=1;iN+1;i+)for(j=0;ji;j+)mi=mi*(varx-xj);mi=t0i*mi;sum=sum+mi;for(i=1;iN+1;i+)printf( +%f*,t0i);for(j=0;ji;j+)printf(x-%f),xj);printf(n);return sum;void main()double varx,fN+1N+1;double xN+1=0.4,0.55,0.8,0.9,1;double yN+1=0.41075,0.57815,0.88811,1.02652,1.17520;checkva

8、ild(x,N);chashang(x,y,f);varx=0.5;if(checkvaild(x,N)=1)chashang(x,y,f);printf(nn牛顿插值结果: P(%f)=%fn,varx,compvalue(f,x,varx);elseprintf(输入的插值节点的x值必须互异!);getch();运行及结果显示：实验三：自适应梯形公式命名（源程序.cpp及工作区.dsw）：autotrap问题：计算的近似值，使得误差（EPS）不超过/autotrap.cpp#include#include#include#define EPS 1e-6double f(double x)r

9、eturn 4/(1+x*x);double AutoTrap(double(*f)(double),double a,double b,double eps)double t2,t1,sum=0.0;int i,k;t1=(b-a)*(f(a)+f(b)/2;printf(T(%4d)=%fn,1,t1);for(k=1;k11;k+)for(i=1;i=pow(2,k-1);i+)sum+=f(a+(2*i-1)*(b-a)/pow(2,k);t2=t1/2+(b-a)*sum/pow(2,k);printf(T(%4d)=%fn,(int)pow(2,k),t2);if(fabs(t2-

10、t1)EPS)break;elset1=t2;sum=0.0;return sum;void main()double s;s=AutoTrap(f,0.0,1.0,EPS);getch();运行及结果显示：实验四：龙贝格算法命名（源程序.cpp及工作区.dsw）：romberg问题：求的近似值，要求误差不超过/romberg.cpp#include #include #include #define N 20#define EPS 1e-6double f(double x)return 4/(1+x*x);double Romberg(double a,double b,double (*

11、f)(double),double eps)double TNN,p,h;int k=1,i,j,m=0; T00=(b-a)/2*(f(a)+f(b);dop=0;h=(b-a)/pow(2,k-1);for(i=1;i=pow(2,k-1);i+) p=p+f(a+(2*i-1)*h/2); T0k=T0k-1/2+p*h/2; for(i=1;i=EPS); k-; while(m=k) for(i=0;i=m;i+) printf(T(%d,%d)=%f ,i,m-i,Tim-i); m+; printf(n); return Tk0;void main()printf(nn积分结果

12、=%f,Romberg(0,1,f,EPS);getch();运行及结果显示：实验五：牛顿切线法问题：求方程在，附近的根（精度）/newton_qxf.cpp#include #include#include #define MaxK 1000#define EPS 0.5e-3double f(double x)return x*x*x-x-1;double f1(double x)return 3*x*x-1;int newton(double (*f)(double), double (*f1)(double), double &x0, double eps)double xk, xk1

13、;int count = 0;printf(ktxkn);printf(-n);xk = x0;printf(%dt%fn, count, xk);doif(*f1)(xk)=0)return 2;xk1 = xk - (*f)(xk) / (*f1)(xk);if (fabs(xk - xk1) eps)count+;xk = xk1;printf(%dt%fn, count, xk);x0 = xk1;return 1;count+;xk = xk1;printf(%dt%fn, count, xk);while(count MaxK);return 0;void main()for(in

14、t i=0;i2;i+)double x=0.6;if(i=1)x=1.5;printf(x0初值为%f:n,x);if (newton(f, f1, x, EPS) = 1)printf(-n);printf(the root is x=%fnnn, x);elseprintf(the method is fail!); getch();实验六：牛顿下山法命名（源程序.cpp及工作区.dsw）：newton_downhill问题：求方程在，附近的根（精度）/newton_downhill.cpp#include #include #include #include #define Et 1e

15、-3/下山因子下界#define E1 1e-3/根的误差限#define E2 1e-3/残量精度double f(double x)return x * x * x - x - 1;double f1(double x)return 3 * x * x - 1;void errormess(int b)char *mess;switch(b)case -1:mess = f(x)的导数为0!;break;case -2:mess = 下山因子已越界，下山处理失败;break;default:mess = 其他类型错误!;printf(the method has faild! becaus

16、e %s, mess);int Newton(double (*f)(double), double (*f1)(double), double &x0)int k = 0;double t;double xk, xk1;double fxk, fxk_temp;printf(k t xk f(xk)n);printf(-n);printf(%-20d, k);xk = x0;printf(%-15f, x0);fxk = (*f)(xk);printf(%-20f, fxk);printf(n);for (k = 1; ; k+)t = 1;while(1)printf(%-10d, k);

17、printf(%-10f, t);if(*f1)(xk) != 0)xk1 = xk - t * (*f)(xk) / (*f1)(xk);elsereturn -1;printf(%-15f, xk1);fxk_temp = (*f)(xk1);printf(%-20f, fxk_temp);if(fabs(fxk_temp) fabs(fxk)t = t / 2;printf(n);if (t Et)return -2;elseprintf(下山成功n);break;if (fabs(xk-xk1)E1)x0 = xk1;return 1;xk = xk1;void main()int b

18、;double x0;x0 = 0.6;b = Newton(f, f1, x0);if (b = 1)printf(nthe root x=%fn, x0);elseerrormess(b);getch();运行及结果显示：实验七：埃特金加速算法命名（源程序.cpp及工作区.dsw）：aitken问题：求方程在，附近的根（精度）/aitken.cpp#include #include #include #define MaxK 100#define EPS 0.5e-3double g(double x)return x * x * x - 1;int aitken (double (*g)

19、(double), double &x, double eps)int k;double xk = x, yk, zk, xk1;printf(k xk yk zk xk+1n);printf(-n);for (k = 0;kMaxK; k+)yk = (*g)(xk);zk = (*g)(yk);xk1 = xk - (yk - xk) * (yk - xk) / (zk - 2 * yk + xk);printf(%-10d%-15f%-15f%-15f%-15fn, k, xk, yk, zk, xk1);if (fabs(yk-xk)=eps)x = xk1;return k + 1;

20、xk = xk1;return -1;void main ()double x = 1.5;int k;k = aitken(g, x, EPS);if (k = -1)printf(迭代次数越界!n);elseprintf(n 经k=%d次迭代，所得方程根为:x=%fn, k, x);getch();运行及结果显示：实验八：正割法问题：求方程在附近的根（精度0.5e-8）/ZhengGe.cpp#include #include #include #define MaxK 1000#define EPS 0.5e-8double f(double x)return x*x*x-x-1;int

21、 ZhengGe(double (*f)(double), double x0, double &x1, double eps)printf(k xk f(xk)n);printf(-n);int k;double xk0, xk, xk1;xk0 = x0;printf(%-10d%-15f%-15fn, 0, x0, (*f)(x0);xk = x1;for (k=1;kMaxK;k+)if(*f)(xk)-(*f)(xk0)=0)return -1;xk1 = xk - (*f)(xk) / ( (*f)(xk) - (*f)(xk0) ) * (xk - xk0);printf(%-1

22、0d%-15f%-15fn, k, xk, (*f)(xk);if (fabs(xk1 - xk)=eps)printf(%-10d%-15f%-15fnn, k + 1, xk1, (*f)(xk1);printf(-nn);x1 = xk1;return 1;xk0 = xk;xk = xk1;return -2;void main ()double x0 = 1, x1 = 2;if (ZhengGe(f, x0, x1, EPS) = 1)printf(the root is x = %fn, x1);elseprintf(the method is fail!);getch();实验

23、九：高斯列选主元算法命名（源程序.cpp及工作区.dsw）：colpivot问题：求解方程组并计算系数矩阵行列式值/colpivot.cpp#include #include #include #define N 3static double aaNN=1,2,-1,1,-1,5,4,1,-2;static double bbN+1=3,0,2;void main()int i,j;double aN+1N+1,bN+1,xN+1,det;double gaussl(double aN+1,double b,double x);for(i=1;i=N;i+)for(j=1;j=N;j+) ai

24、j=aai-1j-1;bi=bbi-1;det=gaussl(a,b,x);if(det!=0)printf(n方程组的解为:);for(i=1;i=N;i+) printf( x%d=%f,i,xi);printf(nn系数矩阵的行列式值%f,det);else printf(nn系数矩阵奇异阵，高斯方法失败 !:);getch();double gaussl(double aN+1,double b,double x)/a传入增广矩阵(有效元素为a1,1.an,n+1),x负责传入迭代初值并返回解向量/(有效值为x1.xn);返回值：系数矩阵行列式值detAdouble det=1.0,F

25、,m,temp;int i,j,k,r; void disp(double aN+1,double b); printf(消元前增广矩阵 :n); disp(a,b); for(k=1;kN;k+) temp=akk; r=k; for(i=k+1;i=N;i+) if(fabs(temp)fabs(aik) temp=aik; r=i;/按列选主元，即确定ik if(ark=0) return 0;/如果aik,k=0,则A为奇异矩阵，停止计算 if(r!=k) for(j=k;j=N;j+) akj+=arj; arj=akj-arj; akj-=arj; bk+=br;br=bk-br;

26、bk-=br;det=-det;/如果ik!=k,则交换A,b第ik行与第k行元素printf(交换第 %d, %d行:n,k,r);disp(a,b); for(i=k+1;i=N;i+) m=aik/akk;for(j=1;j=k;j+)aij=0.0;for(j=k+1;j0;i-) F=0.0; for(j=i+1;j=N;j+) F+=aij*xj;xi=(bi-F)/aii;det*=aNN; return det;void disp(double aN+1,double b)/显示选主元及消元运算中间增广矩阵结果int i,j; for(i=1;i=N;i+) for(j=1;j

27、=N;j+) printf(%10ft,aij); printf(%10fn,bi); 运行及结果显示：实验十：高斯全主元消去算法命名（源程序.cpp及工作区.dsw）：fullpivot问题：利用全主元消去法求解方程组/fullpivot.cpp#include #include #include #define N 3double xN+1;static double aaNN=0.01,2,-0.5,-1,-0.5,2,5,-4,0.5;static double bbN=-5,5,9;void main()int i,j;double aN+1N+1,bN+1,xN+1,det;int

28、 tN+1;/引入列交换保存x向量的各分量的位置顺序double gaussl(double aN+1,double b,double x,int t);for(i=1;i=N;i+)for(j=1;j=N;j+) aij=aai-1j-1;bi=bbi-1;for(i=1;i=1;i-)printf( x%d=%f,ti,xi);if(i1)printf( -);printf(nn系数矩阵的行列式值%f,det);else printf(nn系数矩阵奇异阵，高斯方法失败 !:);getch();double gaussl(double aN+1,double b,double x,int t

29、)/a传入增广矩阵(有效元素为a1,1.an,n+1),x负责传入迭代初值并返回解向量/(有效值为x1.xn);返回值：系数矩阵行列式值detAdouble det=1.0,F,m,temp; int i,j,k,r,s; void disp(double aN+1,double b,int x);/ printf(消元前增广矩阵 :n);/disp(a,b,t);for(k=1;kN;k+) temp=akk;r=k;s=k;for(i=k;i=N;i+) for(j=k;j=N;j+)if(fabs(temp)fabs(aij) temp=aij;r=i;s=j; /选主元，选取ik,jk

30、 if(ars=0) return 0;if(r!=k) for(j=k;j=N;j+) akj+=arj;arj=akj-arj;akj-=arj; bk+=br; br=bk-br; bk-=br; det=-det;printf(交换第 %d, %d行:n,k,r);disp(a,b,t); if(s!=k)for(i=0;i=N;i+) aik+=ais; ais=aik-ais; aik-=ais; tk+=ts; ts=tk=ts; tk-=ts; det=-det; printf(交换第 %d, %d列:n,k,s);disp(a,b,t); for(i=k+1;i=N;i+)

31、m=aik/akk; for(j=1;j=k;j+)aij=0.0; for(j=k+1;j0;i-) F=0.0;for(j=i+1;j=N;j+) F+=aij*xj; xi=(bi-F)/aii; /回代计算 det*=aNN; return det;void disp(double aN+1,double b,int x)/显示选主元及消元运算中间增广矩阵结果int i,j;for(i=1;i=N;i+)for(j=1;j=N;j+)printf(%ft,aij);printf(| x%d = %fn,xi,bi); 运行及结果显示：实验十一：LU分解算法命名（源程序.cpp及工作区.

32、dsw）：LU问题：利用LU法求解方程组/LU.cpp#include#include#include#define N 3static aaNN=3,1,6,2,1,3,1,1,1;static bbN=2,7,4;void main()int i,j;double aN+1N+1,bN+1;void solve(double aN+1,double bN+1);int LU(double aN+1);for(i=1;i=N;i+)for(j=1;j=N;j+)aij=aai-1j-1;bi=bbi-1;if(LU(a)=1)printf(矩阵L如下n);for(i=1;i=N;i+)for

33、(j=1;j=i;j+)if(i=j)printf(%12d ,1);else printf(%12f,aij);printf(n);printf(n矩阵U如下n);for(i=1;i=N;i+)for(j=1;j=N;j+) if(i=j)printf(%12f,aij);else printf(%12s,);printf(n);solve(a,b);for(i=1;i=N;i+)printf(x%d=%f ,i,bi);printf(n);else printf(nthe LU method failed!n); getch();int LU(double aN+1)/对N阶矩A阵进行LU

34、分解，结果L、U存放在A的相应位置 int i,j,k,s; double m,n; for(i=2;i=N;i+) ai1/=a11; for(k=2;k=N;k+) for(j=k;j=N;j+)m=0;for(s=1;sk;s+)m+=aks*asj;akj-=m; if(akk=0) printf(a%d%d=%d ,k,k,0); return -1;/存在元素akk为0 for(i=k+1;i=N;i+)n=0;for(s=1;sk;s+)n+=ais*ask;aik=(aik-n)/akk; return 1;/正常结束void solve(double aN+1,double

35、bN+1)/利用分解的LU求x/回代求解，L和U元素均在A矩阵相应位置；b存放常数列，返回时存放方程组的解double yN+1,F; int i,j; y1=b1; for(i=2;i=N;i+) F=0.0; for(j=1;j0;i-) F=0.0; for(j=N;ji;j-) F+=aij*bj; bi=(yi-F)/aii; 运行及结果显示：实验十二：Guass-Sediel算法命名（源程序.cpp及工作区.dsw）：GS问题：利用GS法求解方程组（精度为）/Guass_sediel.cpp#include iostream.h#include math.h#include std

36、io.h#include#define N 3 /方程的阶数#define MaxK 100 /最大迭代次数#define EPS 0.5e-4 /精度控制static double aaNN=10,-1,-2,-1,10,-2,-1,-1,5;static double bbN=7.2,8.3,4.2;void main()int i,j; double xN+1; double aN+1N+1,bN+1;int GaussSeidel(double aN+1,double b,double x);for(i=1;i=N;i+) for(j=1;j=N;j+) aij=aai-1j-1; bi=bbi-1; xi=0; if(GaussSeidel(a,b,x)=1) printf(the roots is:); for(i=1;i=N;i+) printf( x%d=%f ,i,xi);printf(n); else prin