包括拉格朗日，牛顿插值，高斯，龙贝格，牛顿迭代，牛顿-科特斯，雅克比，秦九昭，幂法，高斯塞德尔 。都是经典的数学算法，希望能开托您的思路。转自kunli.info
1.拉格朗日插值多项式 ，用于离散数据的拟合
- C/C++ code
-
1 #include <stdio.h> 2 #include <conio.h> 3 #include <alloc.h> 4 float lagrange(float *x, float *y, float xx, int n) /*拉格朗日插值算法*/ 5 { 6 int i, j; 7 float *a, yy = 0.0; 8 /*a作为临时变量，记录拉格朗日插值多项式*/ 9 a = (float *) malloc(n * sizeof(float)); 10 for (i = 0; i <= n - 1; i++) { 11 a[i] = y[i]; 12 for (j = 0; j <= n - 1; j++) 13 if (j != i) 14 a[i] *= (xx - x[j]) / (x[i] - x[j]); 15 yy += a[i]; 16 } 17 free(a); 18 return yy; 19 } 20 main() { 21 int i,n; 22 float x[20],y[20],xx,yy; 23 printf("Input n:"); 24 scanf("%d",&n); 25 if(n>=20) { 26 printf("Error!The value of n must in (0,20)."); 27 getch();return 1; 28 } 29 if(n<=0) { 30 printf("Error! The value of n must in (0,20)."); 31 getch(); 32 return 1; 33 }for(i=0;i<=n-1;i++) { 34 printf("x[%d]:",i); 35 scanf("%f",&x[i]); 36 } 37 printf(" "); 38 for(i=0;i<=n-1;i++) { 39 printf("y[%d]:",i); 40 scanf("%f",&y[i]); 41 } 42 printf(" "); 43 printf("Input xx:"); 44 scanf("%f",&xx); 45 yy=lagrange(x,y,xx,n); 46 printf("x=%f,y=%f ",xx,yy); 47 getch(); 48 }
2.牛顿插值多项式，用于离散数据的拟合
- C/C++ code
-
1 #include <stdio.h> 2 #include <conio.h> 3 #include <alloc.h> 4 void difference(float *x, float *y, int n) { 5 float *f; 6 int k, i; 7 f = (float *) malloc(n * sizeof(float)); 8 for (k = 1; k <= n; k++) { 9 f[0] = y[k]; 10 for (i = 0; i < k; i++) 11 f[i + 1] = (f[i] - y[i]) / (x[k] - x[i]); 12 y[k] = f[k]; 13 } 14 return; 15 } 16 main() { 17 int i,n; 18 float x[20],y[20],xx,yy; 19 printf("Input n:"); 20 scanf("%d",&n); 21 if(n>=20) { 22 printf("Error! The value of n must in (0,20)."); 23 getch(); return 1; 24 } 25 if(n<=0) { 26 printf("Error! The value of n must in (0,20)."); 27 getch(); 28 return 1; 29 } 30 for(i=0;i<=n-1;i++) { 31 printf("x[%d]:",i); 32 scanf("%f",&x[i]);}printf(" "); 33 for(i=0;i<=n-1;i++) { 34 printf("y[%d]:",i); 35 scanf("%f",&y[i]); 36 } 37 printf(" "); 38 difference(x,(float *)y,n); 39 printf("Input xx:"); 40 scanf("%f",&xx); 41 yy=y[20]; 42 for(i=n-1;i>=0;i--) yy=yy*(xx-x[i])+y[i]; 43 printf("NewtonInter(%f)=%f",xx,yy); 44 getch();}
3.高斯列主元消去法,求解其次线性方程组
- C/C++ code
-
1 #include<stdio.h> 2 #include <math.h> 3 #define N 20 4 int main() { 5 int n, i, j, k; 6 int mi, tmp, mx; 7 float a[N][N], b[N], x[N]; 8 printf(" Input n:"); 9 scanf("%d", &n); 10 if (n > N) { 11 printf("The input n should in(0,N)! "); 12 getch(); 13 return 1; 14 } 15 if (n <= 0) { 16 printf("The input n should in(0,N)! "); 17 getch(); 18 return 1; 19 } 20 printf("Now input a(i,j),i,j=0...%d: ", n - 1); 21 for (i = 0; i < n; i++) { 22 for (j = 0; j < n; j++) 23 scanf("%f", &a[i][j]); 24 } 25 printf("Now input b(i),i,j=0...%d: ", n - 1); 26 for (i = 0; i < n; i++) 27 scanf("%f", &b[i]); 28 for (i = 0; i < n - 2; i++) { 29 for (j = i + 1, mi = i, mx = fabs(a[i][j]); j < n - 1; j++) 30 if (fabs(a[j][i]) > mx) { 31 mi = j; 32 mx = fabs(a[j][i]); 33 } 34 if (i < mi) { 35 tmp = b[i]; 36 b[i] = b[mi]; 37 b[mi] = tmp; 38 for (j = i; j < n; j++) { 39 tmp = a[i][j]; 40 a[i][j] = a[mi][j]; 41 a[mi][j] = tmp; 42 } 43 } 44 for (j = i + 1; j < n; j++) { 45 tmp = -a[j][i] / a[i][i]; 46 b[j] += b[i] * tmp; 47 for (k = i; k < n; k++) 48 a[j][k] += a[i][k] * tmp; 49 } 50 } 51 x[n - 1] = b[n - 1] / a[n - 1][n - 1]; 52 for (i = n - 2; i >= 0; i--) { 53 x[i] = b[i]; 54 for (j = i + 1; j < n; j++) 55 x[i] -= a[i][j] * x[j]; 56 x[i] /= a[i][i]; 57 } 58 for (i = 0; i < n; i++) 59 printf("Answer: x[%d]=%f ", i, x[i]); 60 getch(); 61 return 0; 62 } 63 #include<math.h> 64 #include<stdio.h> 65 #define NUMBER 20 66 #define Esc 0x1b 67 #define Enter 0x0d 68 float A[NUMBER][NUMBER + 1], ark; int flag, n; 69 exchange(int r,int k); float max(int k); 70 message(); main() {float x[NUMBER]; int r,k,i,j; char celect; clrscr(); printf(" Use Gauss."); printf(" 1.Jie please press Enter."); printf(" 2.Exit press Esc."); celect=getch(); if(celect==Esc) exit(0); printf(" input n="); scanf("%d",&n); printf(" Input matrix A and B:"); for(i=1;i<=n;i++) {printf(" Input a%d1--a%d%d and b%d:",i,i,n,i); for(j=1;j<=n+1;j++) scanf("%f",&A[i][j]);}for(k=1;k<=n-1;k++) {ark=max(k); if(ark==0) {printf(" It'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++) A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];}x[n]=A[n][n+1]/A[n][n]; for( k=n-1;k>=1;k--) {float me=0; for(j=k+1;j<=n;j++) {me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++) {printf(" x%d=%f",i,x[i]);}message();}exchange(int r,int k) {int i; for(i=1;i<=n+1;i++) A[0][i]=A[r][i]; for(i=1;i<=n+1;i++) A[r][i]=A[k][i]; for(i=1;i<=n+1;i++) A[k][i]=A[0][i];}float max( 71 int k) { 72 int i; 73 float temp = 0; 74 for (i = k; i <= n; i++) 75 if (fabs(A[i][k]) > temp) { 76 temp = fabs(A[i][k]); 77 flag = i; 78 } 79 return temp; 80 } 81 message() {printf(" Go on Enter ,Exit press Esc!"); 82 switch(getch()) {case Enter: main(); 83 case Esc: exit(0); 84 default: {printf(" Input error!");message(); 85 } 86 } 87 }
4.龙贝格求积公式，求解定积分
- C/C++ code
-
1 #include<stdio.h> 2 #include<math.h> 3 #define f(x) (sin(x)/x) 4 #define N 20 5 #define MAX 20 6 #define a 2 7 #define b 4 8 #define e 0.00001 9 float LBG(float p, float q, int n) { 10 int i; 11 float sum = 0, h = (q - p) / n; 12 for (i = 1; i < n; i++) 13 sum += f(p+i*h); 14 sum += (f(p) + f(q)) / 2; 15 return (h * sum); 16 } 17 void main() { 18 int i; 19 int n = N, m = 0; 20 float T[MAX + 1][2];T 21 [0][1] = LBG(a, b, n); 22 n *= 2; 23 for (m = 1; m < MAX; m++) { 24 for (i = 0; i < m; i++) 25 T[i][0] = T[i][1]; 26 T[0][1] = LBG(a, b, n); 27 n *= 2; 28 for (i = 1; i <= m; i++) 29 T[i][1] = T[i - 1][1] 30 + (T[i - 1][1] - T[i - 1][0]) / (pow(2, 2 * m) - 1); 31 if ((T[m - 1][1] < T[m][1] + e) && (T[m - 1][1] > T[m][1] - e)) { 32 printf("Answer=%f ", T[m][1]); 33 getch(); 34 return; 35 } 36 } 37 }
-
5.牛顿迭代公式，求方程的近似解
- C/C++ code
-
1 #include<stdio.h> 2 #include<math.h> 3 #include<conio.h> 4 #define N 100 5 #define PS 1e-5 6 #define TA 1e-5 7 float Newton(float(*f)(float), float(*f1)(float), float x0) { 8 float x1, d = 0; 9 int k = 0; 10 do { 11 x1 = x0 - f(x0) / f1(x0); 12 if ((k++ > N) || (fabs(f1(x1)) < PS)) { 13 printf(" Failed!"); 14 getch(); 15 exit(); 16 } 17 d = (fabs(x1) < 1 ? x1 - x0 : (x1 - x0) / x1); 18 x0 = x1; 19 printf("x(%d)=%f ", k, x0); 20 } while ((fabs(d)) > PS && fabs(f(x1)) > TA); 21 return x1; 22 } 23 float f(float x) { 24 return x * x * x + x * x - 3 * x - 3; 25 } 26 float f1(float x) { 27 return 3.0 * x * x + 2 * x - 3; 28 } 29 void main() { 30 float f(float); 31 float f1(float); 32 float x0, y0; 33 printf("Input x0: "); 34 scanf("%f", &x0); 35 printf("x(0)=%f ", x0); 36 y0 = Newton(f, f1, x0); 37 printf(" The root is x=%f ", y0); 38 getch(); 39 }
- 6. 牛顿-科特斯求积公式，求定积分
- C/C++ code
-
1 #include<stdio.h> 2 #include<math.h> 3 int NC(a,h,n,r,f) float (*a)[]; float h; 4 int n, f; 5 float *r; 6 { 7 int nn,i; 8 float ds; 9 if(n>1000||n<2) { 10 if (f) printf(" Faild! Check if 1<n<1000! ",n); 11 return(-1); 12 } 13 if(n==2) {*r=0.5*((*a)[0]+(*a)[1])*(h); 14 return(0); 15 } 16 if (n-4==0) {*r=0; *r=*r+0.375*(h)*((*a)[n-4]+3*(*a)[n-3]+3*(*a)[n-2]+(*a)[n-1]); 17 return(0); 18 } 19 if(n/2-(n-1)/2<=0) nn=n; 20 else nn=n-3; ds=(*a)[0]-(*a)[nn-1]; 21 for(i=2;i<=nn;i=i+2) ds=ds+4*(*a)[i-1]+2*(*a)[i]; *r=ds*(h)/3; 22 if(n>nn) *r=*r+0.375*(h)*((*a)[n-4]+3*(*a)[n-3]+3*(*a)[n-2]+(*a)[n-1]); 23 return(0);} 24 main() { 25 float h,r; 26 int n,ntf,f; 27 int i; 28 float a[16]; 29 printf("Input the x[i](16): "); 30 for(i=0;i<=15;i++) scanf("%d",&a[i]); 31 h=0.2; 32 f=0; 33 ntf=NC(a,h,n,&r,f); 34 if(ntf==0) 35 printf(" R=%f ",r); 36 else 37 printf(" Wrong!Return code=%d ",ntf); 38 getch(); 39 }
7.雅克比迭代，求解方程近似解
- C/C++ code
-
1 #include <stdio.h> 2 #include <math.h> 3 #define N 20 4 #define MAX 100 5 #define e 0.00001 6 int main() { 7 int n; 8 int i, j, k; 9 float t; 10 float a[N][N], b[N][N], c[N], g[N], x[N], h[N]; 11 printf(" Input dim of n:"); 12 scanf("%d", &n); 13 if (n > N) { 14 printf("Faild! Check if 0<n<N! "); 15 getch(); 16 return 1; 17 } 18 if (n <= 0) { 19 printf("Faild! Check if 0<n<N! "); 20 getch(); 21 return 1; 22 } 23 printf("Input a[i,j],i,j=0…%d: ", n - 1); 24 for (i = 0; i < n; i++) 25 for (j = 0; j < n; j++) 26 scanf("%f", &a[i][j]); 27 printf("Input c[i],i=0…%d: ", n - 1); 28 for (i = 0; i < n; i++) 29 scanf("%f", &c[i]); 30 for (i = 0; i < n; i++) 31 for (j = 0; j < n; j++) { 32 b[i][j] = -a[i][j] / a[i][i]; 33 g[i] = c[i] / a[i][i]; 34 } 35 for (i = 0; i < MAX; i++) { 36 for (j = 0; j < n; j++) 37 h[j] = g[j]; 38 { 39 for (k = 0; k < n; k++) { 40 if (j == k) 41 continue; 42 h[j] += b[j][k] * x[k]; 43 } 44 } 45 t = 0; 46 for (j = 0; j < n; j++) 47 if (t < fabs(h[j] - x[j])) 48 t = fabs(h[j] - x[j]); 49 for (j = 0; j < n; j++) 50 x[j] = h[j]; 51 if (t < e) { 52 printf("x_i= "); 53 for (i = 0; i < n; i++) 54 printf("x[%d]=%f ", i, x[i]); 55 getch(); 56 return 0; 57 } 58 printf("after %d repeat , return ", MAX); 59 getch(); 60 return 1; 61 } 62 getch(); 63 }
8.秦九昭算法
- C/C++ code
-
1 #include <math.h> 2 float qin(float a[], int n, float x) { 3 float r = 0; 4 int i; 5 for (i = n; i >= 0; i--) 6 r = r * x + a[i]; 7 return r; 8 } 9 main() { 10 float a[50],x,r=0; 11 int n,i; 12 do {printf("Input frequency:"); 13 scanf("%d",&n); 14 } 15 while(n<1); 16 printf("Input value:"); 17 for(i=0;i<=n;i++) 18 scanf("%f",&a[i]); 19 printf("Input frequency:"); 20 scanf("%f",&x); 21 r=qin(a,n,x); 22 printf("Answer:%f",r); 23 getch(); 24 }
9.幂法
- C/C++ code
-
1 #include<stdio.h> 2 #include<math.h> 3 #define N 100 4 #define e 0.00001 5 #define n 3 6 float x[n] = { 0, 0, 1 }; 7 float a[n][n] = { { 2, 3, 2 }, { 10, 3, 4 }, { 3, 6, 1 } }; 8 float y[n]; 9 main() { 10 int i,j,k; 11 float xm,oxm; 12 oxm=0; 13 for(k=0;k<N;k++) { 14 for(j=0;j<n;j++) { 15 y[j]=0; 16 for(i=0;i<n;i++) y[j]+=a[j][i]*x[i];}xm=0; 17 for(j=0;j<n;j++) 18 if(fabs(y[j])>xm) xm=fabs(y[j]); 19 for(j=0;j<n;j++) y[j]/=xm; 20 for(j=0;j<n;j++) x[j]=y[j]; 21 if(fabs(xm-oxm)<e) { 22 printf("max:%f ",xm); 23 printf("v[i]: "); 24 for(k=0;k<n;k++) 25 printf("%f ",y[k]); 26 break; 27 } 28 oxm=xm; 29 } 30 getch(); 31 }
10.高斯塞德尔
- C/C++ code
-
1 #include<math.h> 2 #include<stdio.h> 3 #define N 20 4 #define M 99 5 float a[N][N]; 6 float b[N]; 7 int main() { 8 int i, j, k, n; 9 float sum, no, d, s, x[N]; 10 printf(" Input dim of n:"); 11 scanf("%d", &n); 12 if (n > N) { 13 printf("Faild! Check if 0<n<N! "); 14 getch(); 15 return 1; 16 } 17 if (n <= 0) { 18 printf("Faild! Check if 0<n<N! "); 19 getch(); 20 return 1; 21 } 22 printf("Input a[i,j],i,j=0…%d: ", n - 1); 23 for (i = 0; i < n; i++) 24 for (j = 0; j < n; j++) 25 scanf("%f", &a[i][j]); 26 printf("Input b[i],i=0…%d: ", n - 1); 27 for (i = 0; i < n; i++) 28 scanf("%f", &b[i]); 29 for (i = 0; i < n; i++) 30 x[i] = 0; 31 k = 0; 32 printf(" k=%dx=", k); 33 for (i = 0; i < n; i++) 34 printf("%12.8f", x[i]); 35 do 36 { 37 k++; 38 if(k>M) 39 { 40 printf (" Error! ”); 41 getch();}break;}no=0.0; 42 for (i = 0; i < n; i++) { 43 s = x[i]; 44 sum = 0.0; 45 for (j = 0; j < n; j++) 46 if (j != i) 47 sum = sum + a[i][j] * x[j]; 48 x[i] = (b[i] - sum) / a[i][i]; 49 d = fabs(x[i] - s); 50 if (no < d) 51 no = d; 52 } 53 printf(" k=%2dx=", k); 54 for (i = 0; i < n; i++) 55 printf("%f", x[i]); 56 } 57 while (no>=0.1e-6); 58 if(no<0.1e-6) { 59 printf(" answer= "); 60 printf(" k=%d",k); 61 for (i=0;i<n;i++) 62 printf(" x[%d]=%12.8f",i,x[i]); 63 }getch(); 64 }