数值分析中的各种公式C代码.docx

《数值分析中的各种公式C代码.docx》由会员分享，可在线阅读，更多相关《数值分析中的各种公式C代码.docx（28页珍藏版）》请在三一办公上搜索。

1、数值分析中的各种公式C 代码二分法 2.2 算法步骤 步骤1： 准备 计算f(x)在有根区间a,b端点处的值f(a)，f(b). 步骤2： 二分 计算f(x)在区间中点 (a+b)/2处的值 f(a+b)/2). 步骤3： 判断 若f(a+b)/2)=0，则(a+b)/2即是根，计算过程结束，否则检验； 若f(a+b)/2)f(a)0，则以(a+b)/2代替b，否则以(a+b)/2代替a. 反复执行步骤2和步骤3，直到区间a,b的长度小于允许误差e，此时中点(a+b)/2即为所 求近似根。 2.3 程序流程图 开始 a,b,e k=0 x0=(a+b)/2 f (a)f(x0)=c N c0

2、Y Y c=0 x0=b 输出x0,k x0=b STOP b-a=e Y 输出x0,k N k+1=k STOP 3 实验结果分析 #includestdio.h void main float a,b,e,x,c; int k=0,n=1; scanf(%f,%f,%f,&a,&b,&e); while(n=1) x=(a+b)/2;c=(x*x*x-x-1)*(a*a*a-a-1); if(c0) b=x; if(b-a=e) printf(%f,%dn,x,k); break; else k+; else if(c=0) printf(%f,%dn,x,k); break; else

3、a=x; if(b-a=e) printf(%f,%dn,x,k); break; else k+; 高斯塞德尔迭代法求方程组解 高斯主元素消去法求方程组解 2.2 算法步骤 高斯塞德尔迭代法： (0)(0)(0)步骤1：给定初值x1,x2,.,xn1 ，精度e，最大容许迭代次数M,令k=1。 nx(1)i=(x(1)(1)i-aijx(j0).aiij=1j0(0)步骤2：对 i=1,2,n 依此计算ei=xi-xixi(1)xi(0)步骤3：求出 e=(0)，若 e，则输出 (i=1,2,.,n )，停止计算。否则执行步 exeiimax1in步骤4：若ke i=1,n,1 (bi-aij

4、xj-j=1i-1j=i+1anijxj)/aijxi,xit y xi-te xi-te i y ee n k 输出x1,x2,.,xn 输出迭代失败标志 结束 高斯主元素消去法： 开始 输入系数矩阵A常数项b及n det-1 k=1,n-1,1 调选列主元子程序 i=k+1,n,1 aikaij-aik/akk j=k+1,n,1 aikaij-aik/akj j bibi-aikbk i detakkdet k 消元过程 bnbn/am i=n-1,1,-1 S-0 j=i+1,n,1 ss=aijbj j bi (bi-s)/aij i detanndet 输出b及det 结束 3 实

5、验结果分析 高斯塞德尔迭代法： 回代过程高斯主元素消去法： 高斯塞德尔迭代法： #include #include/控制输入的格式，如下边setw(6),1.234567控制了小数点后六位 #include #include #define M 10 int n; void Gauss_seidel(double x,double x0,double aM,double b,double eps,int Nmax) int i,j,s=0; double max; while(sNmax) for(i=0;in;i+) xi=bi; for(j=0;jn;j+) if(j!=i) xi=xi-

6、aij*xj; xi=xi/aii; max=fabs(x0-x00); for(i=1;imax) max=fabs(xi-x0i); if(maxeps) break; for(i=0;i=Nmax) cout迭代发散！nendl; else cout原方程组的解为：nendl; for(i=0;in;i+) printf(X%d=%fn,i+1,xi); void main double a1010=5,2,1,-1,4,2,2,-3,10; double b10=-12,20,3; double x10,x0110; double eps;/ 精度 int i,j,s=0,Nmax;

7、cout方程组系数矩阵：endl; n=3; for(i=0;in;i+) for(j=0;jn;j+) coutsetw(6)aij;/同上 coutsetw(6)biendl; cout输入初始向量:endl; for(i=0;ixi; x01i=xi; coutn输入所需精度：eps; cout输入最大迭代次数：Nmax; Gauss_seidel(x,x01,a,b,eps,Nmax); coutendl; 高斯主元素消去法： #include #include #include /在列向量中寻找绝对值最大的项，并返回该项的标号 int FindMax(int p,int N,doub

8、le *A) int i=0,j=0; double max=0.0; for(i=p;imax) j=i; max=fabs(Ai*(N+1)+p); return j; /交换矩阵中的两行 void ExchangeRow(int p,int j,double *A,int N) int i=0; double C=0.0; for(i=0;iN+1;i+) C=Ap*(N+1)+i; Ap*(N+1)+i=Aj*(N+1)+i; Aj*(N+1)+i=C; /上三角变换，A为增广矩阵的指针，N为矩阵的行数。 void uptrbk(double *A,int N) int p=0,k=0

9、,q=0,j=0; double m=0.0; for(p=0;pN-1;p+) /找出该列最大项的标号 j=FindMax(p,N,A); /交换p行和j行 ExchangeRow(p,j,A,N); if(Ap*(N+1)+p=0) printf(矩阵是一个奇异矩阵,没有唯一解!); break; /消去P元素以下的p列内容。 for(k=p+1;kN;k+) m=Ak*(N+1)+p/Ap*(N+1)+p; for(q=p;qN+1;q+) Ak*(N+1)+q=Ak*(N+1)+q-m*Ap*(N+1)+q; printf(n增广矩阵高斯列主元消去后的矩阵为:n); for(j=0;j

10、=0;k-) temp=0.0; for(i=k+1;i256|N=0) printf(输入的数字不再范围之内!); printf(n); return 0; else A=(double*)calloc(N*(N+1),sizeof(double); printf(请输入待求解方程组的增广矩阵(%d行%d列):n,N,N+1); for(i=0;iN*(N+1);i+) scanf(%lf,&Ai); system(cls); printf(方程的增广矩阵为:n); for(i=0;iN*(N+1);i+) if(i%(N+1)=0) printf(n); printf(%lft,Ai);

11、uptrbk(A,N); /上三角变换 X=backsub(A,N); /回代函数 printf(nn方程组的解为:n); for(i=0;iN;i+) printf(X(%d)= %lfn,i+1,Xi); free(A); free(X); exit(0); 二次或者三次样条插值 #include math.h #include stdio.h #include #include iostream.h #include GaussRemove.h #define BOOL int #define FALSE 0 #define TRUE 1 /计算分段线性插值的系数 void Linear

12、Coe(double * pd_X, double * pd_Y,double * pd_Coe, int i) /计算分段二次多项式插值的系数 void QuadrateCoe(double * pd_X, double * pd_Y,double * pd_Coe, int i) for(int j = 0; j n+1; j+) GaussRemove(pd_MatrixA,pd_VectorB,pd_Coe,n+1); pd_MatrixAj*3 = pd_Xi+j-1*pd_Xi+j-1; pd_MatrixAj*3 + 1 = pd_Xi+j-1; pd_MatrixAj*3 +

13、2 = 1; pd_VectorB0 = pd_Yi-1; pd_VectorB1 = pd_Yi; pd_VectorB2 = pd_Yi+1; int n = 2; double * pd_MatrixA = new double(n+1)*(n+1); double * pd_VectorB = new doublen+1; pd_Coe0 = (-pd_Xi*pd_Yi+1+pd_Xi+1*pd_Yi)/(pd_Xi+1-pd_Xi); pd_Coe1 = (pd_Yi+1-pd_Yi)/(pd_Xi+1-pd_Xi); /计算分段三次多项式插值的系数 void ThriceCoe(d

14、ouble * pd_X, double * pd_Y,double * pd_Coe, int i) /计算三次样条插值函数的系数组M void SplineCoe(double * pd_X, double * pd_Y, int N, double * pd_h, double * pd_M) int n = N -1; double * pd_Alpha = new doublen; /系数数组Alpha double * pd_Beta = new doublen; /系数数组Beta double * pd_Gama = new doublen; /系数数组Gama int i,j

15、; int n = 3; double * pd_MatrixA = new double(n+1)*(n+1); double * pd_VectorB = new doublen+1; pd_VectorB0 = pd_Yi-2; pd_VectorB1 = pd_Yi-1; pd_VectorB2 = pd_Yi; pd_VectorB3 = pd_Yi+1; for(int j = 0; j n+1; j+) GaussRemove(pd_MatrixA,pd_VectorB,pd_Coe,n+1); pd_MatrixAj*4 = pd_Xi+j-2*pd_Xi+j-2*pd_Xi+

16、j-2; pd_MatrixAj*4 + 1 = pd_Xi+j-2*pd_Xi+j-2; pd_MatrixAj*4 + 2 = pd_Xi+j-2; pd_MatrixAj*4 + 3 = 1; for(i = 0; i n; i+) pd_Alphai = pd_h0/(pd_h0 + pd_hi); pd_Gamai = 1 - pd_Alphai; pd_Betai = 6 * (pd_Y1-pd_Y0)/pd_h0) - (pd_Yi+1-pd_Yi)/pd_hi) / /计算Alpha,Beta,Gama (pd_h0 + pd_hi); double * pd_MatrixA

17、= new doublen*n; for(i = 0; i n; i+) for(j = 0; j n; j+) if (i = j) else if (j = (i+1) else if (j = (i-1) else pd_MatrixAi*n+j = 0; pd_MatrixAi*n+j = pd_Gamai; pd_MatrixAi*n+j = pd_Alphai; pd_MatrixAi*n+j = 2; pd_MatrixAn-1 = pd_Gama0; pd_MatrixA(n-1)*n = pd_Alphan-1; / ShowMatrix(pd_MatrixA,n,n); /

18、 ShowVector(pd_Beta,n); delete pd_Alpha; delete pd_Beta; delete pd_Gama; delete pd_MatrixA; GaussRemove(pd_MatrixA,pd_Beta,&(pd_M1),N-1); pd_M0 = pd_MN-1; 复化梯形公式和复化辛普森公式 对应程序 实验一的源程序： #include #include double fun1(double x,double y,int n) int i; double sum=0,m; for(i=1;i=n-1;i+) /求f(x(k)的和 m=x+i*y;

19、sum=sum+sqrt(1+pow(2.718,m); return sum; double get(double x) return sqrt(1+pow(2.718,x); double fun2(double x,double y,int n) int i; double sum1,sum2=0,ff; sum1=get(x+y/2); for(i=1;i=n-1;i+)/求f(x(k)的和 sum1+=get(x+i*y+y/2); sum2+=get(x+i*y); ff=4*sum1+2*sum2; return ff; void main double a=0,b=2;/给定的

20、区间 int n; while(1) cout0):; cinn; double h; double I1,f10,f11,f12; double I2,f20,f21,f22; h=(b-a)/n;/步长 f10=get(a); f11=fun1(a,h,n); f12=get(b); f20=f10; f21=fun2(a,h,n); f22=f12; I1=h/2*(f10+2*f11+f12);/复化梯形公式 I2=h/6*(f20+f21+f22); /复化辛浦生公式 cout等分成n个空间endl; cout复化梯形公式计算结果I=I1endl; cout复化辛浦生公式计算结果I=

21、I2endl; 实验内容二的源程序： #include #include double get(double x) if(x=0) return 1; else return sin(x)/x; double fun(double x,double y,int n) int i; double sum1,sum2=0,ff; sum1=get(x+y/2); for(i=1;i=n-1;i+)/求f(x(k)的和 sum1+=get(x+i*y+y/2); sum2+=get(x+i*y); ff=4*sum1+2*sum2; return ff; void main double a=0,b

22、=1;/给定的区间 int n;/等分成n个空间 while(true) cout0):; cinn; if(n=-1) break; double h; double I1,f11,f12,f13; h=(b-a)/n;/步长 f11=get(a); f12=fun(a,h,n); f13=get(b); I1=h/6*(f11+f12+f13);/复化辛浦生公式 cout等分成n个区间:endl; cout复化辛浦生公式计算结果:I=I1endl; 10个重要的算法C语言实现源代码：拉格朗日，牛顿插值，高斯，龙贝格，牛顿迭代，牛顿-科特斯，雅克比，秦九昭，幂法，高斯塞德尔 (转) 1.拉格

23、朗日插值多项式 ，用于离散数据的拟合 C/C+ code#include #include #include float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法*/ int i,j; float *a,yy=0.0; /*a作为临时变量，记录拉格朗日插值多项式*/ a=(float *)malloc(n*sizeof(float); for(i=0;i=n-1;i+) ai=yi; for(j=0;j=20) printf(Error!The value of n must in (0,20).); getch;return

24、1; if(n=0) printf(Error! The value of n must in (0,20).); getch; return 1; for(i=0;i=n-1;i+) printf(x%d:,i); scanf(%f,&xi); printf(n); for(i=0;i=n-1;i+) printf(y%d:,i);scanf(%f,&yi); printf(n); printf(Input xx:); scanf(%f,&xx); yy=lagrange(x,y,xx,n); printf(x=%f,y=%fn,xx,yy); getch; 2.牛顿插值多项式，用于离散数据

25、的拟合 C/C+ code#include #include #include void difference(float *x,float *y,int n) float *f; int k,i; f=(float *)malloc(n*sizeof(float); for(k=1;k=n;k+) f0=yk; for(i=0;i=20) printf(Error! The value of n must in (0,20).); getch; return 1; if(n=0) printf(Error! The value of n must in (0,20).);getch; ret

26、urn 1; for(i=0;i=n-1;i+) printf(x%d:,i); scanf(%f,&xi); printf(n); for(i=0;i=0;i-) yy=yy*(xx-xi)+yi; printf(NewtonInter(%f)=%f,xx,yy); getch; 3.高斯列主元消去法,求解其次线性方程组 C/C+ code#include #include #define N 20 int main int n,i,j,k; int mi,tmp,mx; float aNN,bN,xN; printf(nInput n:); scanf(%d,&n); if(nN) pri

27、ntf(The input n should in(0,N)!n); getch; return 1; if(n=0) printf(The input n should in(0,N)!n); getch; return 1; printf(Now input a(i,j),i,j=0.%d:n,n-1); for(i=0;in;i+) for(j=0;jn;j+) scanf(%f,&aij); printf(Now input b(i),i,j=0.%d:n,n-1); for(i=0;in;i+) scanf(%f,&bi); for(i=0;in-2;i+) for(j=i+1,mi

28、=i,mx=fabs(aij);jmx) mi=j; mx=fabs(aji); if(imi) tmp=bi;bi=bmi;bmi=tmp; for(j=i;jn;j+) tmp=aij; aij=amij; amij=tmp; for(j=i+1;jn;j+) tmp=-aji/aii; bj+=bi*tmp; for(k=i;k=0;i-) xi=bi; for(j=i+1;jn;j+) xi-=aij*xj; xi/=aii; for(i=0;in;i+) printf(Answer:n x%d=%fn,i,xi); getch; return 0; #include #include

29、 #define NUMBER 20 #define Esc 0x1b #define Enter 0x0d float ANUMBERNUMBER+1 ,ark; int flag,n; exchange(int r,int k); float max(int k); message; main float xNUMBER; int r,k,i,j; char celect; clrscr; printf(nnUse Gauss.); printf(nn1.Jie please press Enter.); printf(nn2.Exit press Esc.); celect=getch;

30、 if(celect=Esc) exit(0); printf(nn input n=); scanf(%d,&n); printf( nnInput matrix A and B:); for(i=1;i=n;i+) printf(nnInput a%d1-a%d%d and b%d:,i,i,n,i); for(j=1;j=n+1;j+) scanf(%f,&Aij); for(k=1;k=n-1;k+) ark=max(k); if(ark=0) printf(nnIts wrong!);message; else if(flag!=k) exchange(flag,k); for(i=

31、k+1;i=n;i+) for(j=k+1;j=1;k-) float me=0; for(j=k+1;j=n;j+) me=me+Akj*xj; xk=(Akn+1-me)/Akk; for(i=1;i=n;i+) printf( nnx%d=%f,i,xi); message; exchange(int r,int k) int i; for(i=1;i=n+1;i+) A0i=Ari; for(i=1;i=n+1;i+) Ari=Aki; for(i=1;i=n+1;i+) Aki=A0i; float max(int k) int i; float temp=0; for(i=k;itemp) temp=fabs(Aik); flag=i; return temp; message printf(nn Go on Enter ,Exit press Esc!); switch(getch) case Enter: main; case Esc: