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

数值分析中的各种公式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 c<0 Y Y c=0 x0=b 输出x0,k x0=b STOP b-a<=e Y 输出x0,k N k+1=k STOP 3 实验结果分析 #include"stdio.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(c<0) 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 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=1j¹0(0)步骤2：对 i=1,2,n 依此计算ei=xi-xixi(1)®xi(0)步骤3：求出 e=(0)，若 e<，则输出 (i=1,2,.,n )，停止计算。否则执行步 exeiimax1£i£n步骤4：若k<M, k+1®k,转步骤2继续迭代。若k³M,表明迭代失败，停止计算。 高斯主元素消去法： 步骤1：det¬1，k=1,2,n-1,从第2步做到第7步； 步骤2： 按列主元素aik,k=maxaik k£i£n步骤3：如果 aik,k=0，则 det¬0，计算停止。 步骤4：若 ik=k，转步骤5，否则换行： akj«aik,j,(j=k,k+1,.n),bk«bik,det¬det 步骤5：计算乘数mik , aik¬mik=aik/akk (i=k_1,n) 步骤6：消元计算 aij¬aij-mikakj,(i,j=k+1,.,n)bi¬bi-mikbi,(i=k+1,.,n)步骤7：det¬akkdet; 步骤8：若 ann=0,det¬0计算停止。否则 det¬anndet bn¬bn/ann步骤9：回代求解 bi¬(aij-j=i+1åab)aijjn (i=n-1,1) ii2.3 程序流程图 高斯塞德尔迭代法： 开始 输入aij,bi,xi及e,m k=1,M,1 0->e i=1,n,1 (bi-åaijxj-j=1i-1j=i+1åanijxj)/aij®xi,xi®t y xi-t£e xi-t®e i y e<e n k 输出x1,x2,.,xn 输出迭代失败标志 结束 高斯主元素消去法： 开始 输入系数矩阵A常数项b及n det<-1 k=1,n-1,1 调选列主元子程序 i=k+1,n,1 aik¬aij-aik/akk j=k+1,n,1 aik¬aij-aik/akj j bi¬bi-aikbk i det¬akkdet k 消元过程 bn¬bn/am i=n-1,1,-1 S<-0 j=i+1,n,1 s¬s=aijbj j bi¬ (bi-s)/aij i det¬anndet 输出b及det 结束 3 实验结果分析 高斯塞德尔迭代法： 回代过程高斯主元素消去法： 高斯塞德尔迭代法： #include<iostream.h> #include<iomanip.h>/控制输入的格式，如下边setw(6),1.234567控制了小数点后六位 #include<stdio.h> #include<math.h> #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(s<Nmax) for(i=0;i<n;i+) xi=bi; for(j=0;j<n;j+) if(j!=i) xi=xi-aij*xj; xi=xi/aii; max=fabs(x0-x00); for(i=1;i<n;i+) if(fabs(xi-x0i)>max) max=fabs(xi-x0i); if(max<eps) break; for(i=0;i<n;i+) x0i=xi; s+; if(s>=Nmax) cout<<"迭代发散！n"<<endl; else cout<<"原方程组的解为：n"<<endl; for(i=0;i<n;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; cout<<"方程组系数矩阵："<<endl; n=3; for(i=0;i<n;i+) for(j=0;j<n;j+) cout<<setw(6)<<aij;/同上 cout<<setw(6)<<bi<<endl; cout<<"输入初始向量:"<<endl; for(i=0;i<n;i+) printf("X%d=",i+1); cin>>xi; x01i=xi; cout<<"n输入所需精度："<<endl; cin>>eps; cout<<"输入最大迭代次数："<<endl; cin>>Nmax; Gauss_seidel(x,x01,a,b,eps,Nmax); cout<<endl; 高斯主元素消去法： #include<stdio.h> #include<stdlib.h> #include<math.h> /在列向量中寻找绝对值最大的项，并返回该项的标号 int FindMax(int p,int N,double *A) int i=0,j=0; double max=0.0; for(i=p;i<N;i+) if(fabs(Ai*(N+1)+p)>max) 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;i<N+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,q=0,j=0; double m=0.0; for(p=0;p<N-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;k<N;k+) m=Ak*(N+1)+p/Ap*(N+1)+p; for(q=p;q<N+1;q+) Ak*(N+1)+q=Ak*(N+1)+q-m*Ap*(N+1)+q; printf("n增广矩阵高斯列主元消去后的矩阵为:n"); for(j=0;j<N*(N+1);j+) if(j%(N+1)=0) printf("n"); printf("%lft",Aj); /下面是回代函数 double* backsub(double *A,int N) double* X=NULL,temp=0.0; int k=0,i=0; X=(double*)malloc(N*sizeof(double); XN-1=A(N-1)*(N+1)+N/A(N-1)*(N+1)+N-1; for(k=N-2;k>=0;k-) temp=0.0; for(i=k+1;i<N;i+) temp=temp+Ak*(N+1)+i*Xi; Xk=(Ak*(N+1)+N-temp)/Ak*(N+1)+k; return X; main int N=0,i=0; double *A=NULL,*X=NULL; printf("n请输入待求解方程组的增广矩阵的行数:"); scanf("%d",&N); if(N>256|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;i<N*(N+1);i+) scanf("%lf",&Ai); system("cls"); printf("方程的增广矩阵为:n"); for(i=0;i<N*(N+1);i+) if(i%(N+1)=0) printf("n"); printf("%lft",Ai); uptrbk(A,N); /上三角变换 X=backsub(A,N); /回代函数 printf("nn方程组的解为:n"); for(i=0;i<N;i+) printf("X(%d)= %lfn",i+1,Xi); free(A); free(X); exit(0); 二次或者三次样条插值 #include "math.h" #include "stdio.h" #include <iomanip.h> #include "iostream.h" #include "GaussRemove.h" #define BOOL int #define FALSE 0 #define TRUE 1 /计算分段线性插值的系数 void LinearCoe(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 + 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(double * 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; 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+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 = 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); / 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<iostream.h> #include<math.h> 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; 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;/给定的区间 int n; while(1) cout<<"请输入要等分的子区间个数(n>0):" cin>>n; 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="<<I1<<endl; cout<<"复化辛浦生公式计算结果I="<<I2<<endl; 实验内容二的源程序： #include<iostream.h> #include<math.h> 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=1;/给定的区间 int n;/等分成n个空间 while(true) cout<<"请输入要等分的子区间个数(n>0):" cin>>n; 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="<<I1<<endl; 10个重要的算法C语言实现源代码：拉格朗日，牛顿插值，高斯，龙贝格，牛顿迭代，牛顿-科特斯，雅克比，秦九昭，幂法，高斯塞德尔 (转) 1.拉格朗日插值多项式 ，用于离散数据的拟合 C/C+ code#include <stdio.h> #include <conio.h> #include <alloc.h> 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<=n-1;j+) if(j!=i) ai*=(xx-xj)/(xi-xj); yy+=ai; free(a); return yy; main int i,n; float x20,y20,xx,yy; printf("Input n:"); scanf("%d",&n); if(n>=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; 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.牛顿插值多项式，用于离散数据的拟合 C/C+ code#include <stdio.h> #include <conio.h> #include <alloc.h> 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<k;i+) fi+1=(fi-yi)/(xk-xi); yk=fk; return; main int i,n; float x20,y20,xx,yy; printf("Input n:"); scanf("%d",&n); if(n>=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; 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"); difference(x,(float *)y,n); printf("Input xx:"); scanf("%f",&xx); yy=y20; for(i=n-1;i>=0;i-) yy=yy*(xx-xi)+yi; printf("NewtonInter(%f)=%f",xx,yy); getch; 3.高斯列主元消去法,求解其次线性方程组 C/C+ code#include<stdio.h> #include <math.h> #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(n>N) printf("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;i<n;i+) for(j=0;j<n;j+) scanf("%f",&aij); printf("Now input b(i),i,j=0.%d:n",n-1); for(i=0;i<n;i+) scanf("%f",&bi); for(i=0;i<n-2;i+) for(j=i+1,mi=i,mx=fabs(aij);j<n-1;j+) if(fabs(aji)>mx) mi=j; mx=fabs(aji); if(i<mi) tmp=bi;bi=bmi;bmi=tmp; for(j=i;j<n;j+) tmp=aij; aij=amij; amij=tmp; for(j=i+1;j<n;j+) tmp=-aji/aii; bj+=bi*tmp; for(k=i;k<n;k+) ajk+=aik*tmp; xn-1=bn-1/an-1n-1; for(i=n-2;i>=0;i-) xi=bi; for(j=i+1;j<n;j+) xi-=aij*xj; xi/=aii; for(i=0;i<n;i+) printf("Answer:n x%d=%fn",i,xi); getch; return 0; #include<math.h> #include<stdio.h> #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; 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("nnIt's wrong!");message; else if(flag!=k) exchange(flag,k); for(i=k+1;i<=n;i+) for(j=k+1;j<=n+1;j+) Aij=Aij-Akj*Aik/Akk; xn=Ann+1/Ann; for( k=n-1;k>=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;i<=n;i+) if(fabs(Aik)>temp) temp=fabs(Aik); flag=i; return temp; message printf("nn Go on Enter ,Exit press Esc!"); switch(getch) case Enter: main; case Esc: