线性方程组求解

高斯消元法

void input(double a[MAXSIZE][MAXSIZE+1],long n);
5 void output(double x[MAXSIZE],long n);
6 void main(void) {
7 double a[MAXSIZE][MAXSIZE+1],x[MAXSIZE],s;
8 long n,i,j,k;
9 printf(”\n请输入原方程组的阶数：”);
10 scanf(”%ld”,&n);
11 input(a,n);
12 for(k=0; k<=n−2; k++)
13 for(i=k+1; i<=n−1; i++) {
14 a[i][k]/=−a[k][k];
15 for(j=k+1; j<=n; j++)
16 a[i][j]+=a[i][k]*a[k][j];
17 }
18 for(k=n−1; k>=0; k−−) {
19 s=0;
20 for(j=k+1; j<=n−1; j++) s+=a[k][j]*x[j];
21 x[k]=(a[k][n]−s)/a[k][k];
22 }
23 output(x,n);
24 }
26 void input(double a[MAXSIZE][MAXSIZE+1],long n) {
27 long i,j;
28 printf(”\n请输入原方程组的增广矩阵：\n”);
29 for(i=1; i<=n; i++)
30 for(j=1; j<=n+1; j++)
31 scanf(”%lf”,&a[i−1][j−1]);
32 }
33 void output(double x[MAXSIZE],long n) {
34 long k;
35 printf(”\n原方程组的解为：\n”);
36 for(k=1; k<=n; k++)
37 printf(” %lf”,x[k−1]);
38 }

列主元高斯消去法

float *colPivot( float *c,int n ) {
7 int i,j,t,k;
8 float *x,p;
9 x=new float[n*sizeof(float)];
10 for( i=0; i<=n−2; i++) {
11 k=i;
12 for(j=i+1; j<=n−1; j++)
13 if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))
14 k=j;
15 if(k!=i)
16 for( j=i; j<=n; j++ ) {
17 p=*(c+i*(n+1)+j);
18 *(c+i*(n+1)+j)=*(c+k*(n+1)+j);
19 *(c+k*(n+1)+j)=p;
20 }
21 for( j=i+1; j<=n−1; j++ ) {
22 p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));
23 for( t=i; t<=n; t++ )
24 *(c+j*(n+1)+t)−=p*(*(c+i*(n+1)+t));
25 }
26 }
27 for( i=n−1; i>=0; i−−) {
28 for( j=n−1; j>=i+1; j−−)
29 (*(c+i*(n+1)+n))−=x[j]*(*(c+i*(n+1)+j));
30 x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));
31 }
32 return x;
33 }
35 int main() {
36 int i;
37 float *x;
38 float c[3][4] = {0.101,2.304,3.555,1.183,
39 −1.347,3.712,4.623,2.137,
40 −2.835,1.072,5.643,3.035
41 };
42 float *ColPivot(float *,int);
43 x=colPivot(c[0],3);
44 for( i=0; i<=2; i++ )
45 cout<<”x(”<<i<<”)=”<<x[i]<<endl;
46 return 0;
47 }

杜利特尔 (Doolittle) 分解

#define MAX_N 20
7 int main(int argc, char*argv[]) {
8 int n; // 未知数个数
9 int i, j, k;
10 static double a[MAX_N][MAX_N], b[MAX_N], x[MAX_N], y[MAX_N];
11 static double l[MAX_N][MAX_N], u[MAX_N][MAX_N];
12 printf(”\nInput n value(dim of Ax=b): ”); //输入系数矩阵维度
13 scanf(”%d”, &n);
14 if(n >MAX_N) {
15 printf(”The input n is larger than MAX_N, please redefine the MAX_N.\n”);
16 return 1;
17 }
18 if(n <= 0) {
19 printf(”Please input a number between 1 and %d.\n”, MAX_N);
20 return 1;
21 }
22 printf(”Now input the matrix a(i, j), i, j = 0, ..., %d:\n”, n−1); //输入 系数矩阵
23 for (i=0; i<n; i++)
24 for (j=0; j<n; j++)
25 scanf(”%lf”, &a[i][j]);
26 printf(”Now input the matrix b(i), i = 0, ..., %d:\n”, n−1); //输入常数项
27 for(i=0; i<n; i++)
28 scanf(”%lf”, &b[i]);
29 for(i=0; i<n; i++)
30 l[i][i] = 1;
31 for(k=0; k<n; k++) {
32 for(j=k; j<n; j++) { // dolittle分解
33 u[k][j]=a[k][j];
34 for(i=0; i<=k−1; i++)
35 u[k][j]−=(l[k][i]*u[i][j]);
36 printf(”%f\n”, u[k][j]);
37 }
38 for(i=k+1; i<n; i++) {
39 l[i][k]=a[i][k];
40 for(j=0; j<=k−1; j++)
41 l[i][k]−=(l[i][j]*u[j][k]);
42 l[i][k]/=u[k][k];
43 printf(”%f\n”, l[i][k]);
44 }
45 }
47 for(i=0; i<n; i++) { // 解Ly = b
48 y[i] = b[i];
49 for(j=0; j<=i−1; j++)
50 y[i] −= (l[i][j]*y[j]);
51 }
53 for(i=n−1; i>=0; i−−){ // 解UX = Y
54 x[i]=y[i];
55 for(j=i+1; j<n; j++)
56 x[i] −= (u[i][j]*x[j]);
57 x[i]/=u[i][i];
58 }
60 printf(”Solve...x_i = \n”); // 输出结果
61 for(i=0; i<n; i++)
62 printf(”%f\n”, x[i]);
63 system(”pause”);
64 return 0;
65 }

克劳特 (Crout) 分解

#define N 20//自定义N=20
5 int main() { //主函数
6 int i,j,k;
7 int size;
8 float a[N][N],l[N][N],u[N][N];
9 float b[N],x[N],y[N];//定义变量
10 printf(”\t\t\tCrout分解法解方程组\n”);
11 printf(”请输入方阵A的n：”);
12 scanf(”%d”,&size);
13 printf(”\n”);
14 printf(”请输入方程组的系数:\n”);
15 for(i=0; i<size; i++) {
16 for(j=0; j<size; j++) {
17 scanf(”%f”,&a[i][j]);//输入方程组系数矩阵a[][]
18 }
19 }
20 printf(”\n请输入方程组的y:\n”);
21 for(i=0; i<size; i++)
22 scanf(”%f”,&b[i]);//输入结果矩阵b[]
23 printf(”\n方阵A[][]为：\n”);
24 for(i=0; i<size; i++) {
25 for(j=0; j<size; j++) {
26 printf(”%f ”,a[i][j]);//输出a[][]
27 }
28 printf(”\n”);
29 }
31 printf(”\n方程组y为:\n”);
32 for(i=0; i<size; i++)
33 printf(”%f ”,b[i]);//输出b[]
34 printf(”\n”);
35 for(i=0; i<size; i++) {
36 u[i][i]=1;//定初始值 令u[i][i]=1
37 }
38 for(i=0; i<size; i++)
39 for(j=i+1; j<size; j++) {
40 l[i][j]=0;//定初始值 令l[i][j]=0
41 }
42 for(j=0; j<size; j++)
43 for(i=j+1; i<size; i++) {
44 u[i][j]=0;//定初始值 令u[i][j]=0
45 }
46 l[0][0]=a[0][0];
47 for(i=1; i<size; i++) {
48 l[i][0]=a[i][0];//计算第一行的l[][]
49 u[0][i]=a[0][i]/l[0][0];//计算第一列的u[][]
50 }
51 for(i=1; i<size−1; i++) {
52 for(j=1; j<=i; j++) { //计算第2行到第size−1行的l[][]
53 l[i][j]=a[i][j];
54 for(k=0; k<j; k++) {
55 l[i][j]=l[i][j]−l[i][k]*u[k][j];
56 }
57 }
58 printf(”\n”);
59 for(j=i+1; j<size; j++) { //计算第2行到第size行的u[][]
60 u[i][j]=a[i][j];
61 for(k=0; k<=i−1; k++) {
62 u[i][j]=u[i][j]−l[i][k]*u[k][j];
64 }
65 u[i][j]=u[i][j]/l[i][i];
66 }
67 printf(”\n”);
68 }
69 for(j=1; j<size; j++) { //计算第size行的l[][]
70 l[size−1][j]=a[size−1][j];
71 for(k=0; k<=j−1; k++) {
72 l[size−1][j]=l[size−1][j]−l[size−1][k]*u[k][j];
73 }
74 }
75 printf(”\n”);
76 printf(”输出矩阵L[i][j]\n”);
77 for(i=0; i<size; i++) {
78 for(j=0; j<size; j++) {
79 printf(”%f”,l[i][j]);
80 printf(” ”);//输出下三角矩阵l[][]
81 }
82 printf(”\n”);
83 }
84 printf(”输出矩阵U[i][j]\n”);
85 for(i=0; i<size; i++) {
86 for(j=0; j<size; j++) {
87 printf(”%f”,u[i][j]);
88 printf(” ”);//输出单位上三角矩阵u[][]
89 }
90 printf(”\n”);
91 }
92 y[0]=b[0]/l[0][0];//给y[0]初始值
93 for(i=1; i<size; i++) { //计算y[i]的值
94 y[i]=b[i];
95 for(k=0; k<=i−1; k++) {
96 y[i]=y[i]−l[i][k]*y[k];//计算公式
97 }
98 y[i]=y[i]/l[i][i];
99 }
100 printf(”\n”);
101 printf(”y值:\n”);
102 for(i=0; i<size; i++)
103 printf(”y[%d]=%f ”,i+1,y[i]);//输出y[i]的结果
104 printf(”\n\n”);
105 printf(”x的值:\n”);
106 x[size−1]=y[size−1];//给x[size−1]赋值
107 for(i=size−2; i>=0; i−−) {
108 x[i]=y[i];
109 for(k=i+1; k<size; k++) {
110 x[i]=x[i]−u[i][k]*x[k];//计算x[]
111 }
112 }
113 for(i=0; i<size; i++) {
114 printf(”x[%d]=%f\n”,i+1,x[i]);//输出x[i]的结果
115 }
116 }

Cholesky 分解

class Cholesky {
12 private:
13 int i,j,k,n;
14 double sum,*b,*d,*x,**a,eps;
15 public:
16 void choleskyInput();
17 void choleskyDecomposition();
18 void choleskyOutput();
19 ~Cholesky() {
20 delete []b;
21 delete []d;
22 delete []x;
23 for(i=0; i<n; i++) {
24 delete [] a[i];
25 }
26 delete []a;
27 }
29 };
31 void main() {
32 Cholesky solution;
33 solution.choleskyInput();
34 solution.choleskyDecomposition();
35 solution.choleskyOutput();
36 }
38 void Cholesky::choleskyInput() {
39 cout<<”输入方程的个数”;
40 cin>>n;
41 b=new double[n];
42 d=new double[n];
43 x=new double[n];
44 a=new double*[n];
46 for(i=0; i<n; i++) {
47 a[i] = new double[n];
48 }
50 for(i=0; i<n; i++)
51 for(j=0; j<n; j++) {
52 cout<<”/n输入a[”<<i<<”][”<<j<<”]=”;
53 cin>>a[i][j];
54 }
55 for(i=0; i<n; i++)
56 for(j=0; j<n; j++) {
57 if(a[i][j] != a[j][i]) {
58 cout<<”/n系数矩阵不对称.失败...”<<endl;
59 exit(0);
60 }
61 }
62 for(i=0; i<n; i++) {
64 cout<<”/n输入 b[”<<i<<”]=”;
65 cin>>b[i];
66 }
67 cout<<”/n输入最小主元素”;
68 cin>>eps; //输入段结束
69 }
70 void Cholesky::choleskyDecomposition() {
71 for(i=0; i<n; i++)
72 for(j=0; j<n; j++) {
73 sum=a[i][j];
74 for(k=0; k<i; k++) {
75 sum −= a[i][k]*a[j][k];
76 }
77 if( i==j) {
78 if(sum <= 0) {
79 cout<<”/n矩阵非正定.失败...”<<endl;
80 exit(0);
81 }
82 d[i]=sqrt(sum);
83 } else {
84 a[j][i] = sum/d[i];
85 }
87 }
89 for(i=0; i<n; i++) {
90 sum=b[i];
91 for(k=0; k<i; k++) {
92 sum −= a[i][k]*x[k];
93 }
94 x[i] = sum/d[i];
95 }
97 for(i=(n−1); i>=0; i−−) {
98 sum = x[i];
99 for(k=(i+1); k<n; k++) {
100 sum −= a[k][i]*x[k];
101 }
102 x[i] = sum/d[i];
103 }
104 }
105 void Cholesky::choleskyOutput() {
106 cout<<”/n:结果是：”<<endl;
107 for(i=0; i<n; i++) {
108 cout<<”x[”<<i<<”]=”<<x[i]<<endl;
109 }
110 }

追赶法

#include <stdio.h>
2 #include <conio.h>
3 #include <dos.h>
5 void main() {
6 char play;
7 int i,n;
8 float a[100],b[100],c[100],d[100];
9 float u[100],l[100],y[100],x[100];
11 do {
12 printf(”Please input the square distance set of piece n:\r\n”);
13 scanf(”%d”,&n);
15 printf(”Please input a2−−an:\n”);
16 for(i=2; i<=n; i++) {
17 printf(”a%d=”,i);
18 scanf(”%f”,&a[i]);
19 }
21 printf(”Please input b1−−bn:\n”);
22 for(i=1; i<=n; i++) {
23 printf(”b%d=”,i);
24 scanf(”%f”,&b[i]);
25 }
27 printf(”Please input c1−−c(n−1):\n”);
28 for(i=1; i<n; i++) {
29 printf(”c%d=”,i);
30 scanf(”%f”,&c[i]);
31 }
33 printf(”Please input d1−−dn:\n”);
34 for(i=1; i<=n; i++) {
35 printf(”d%d=”,i);
36 scanf(”%f”,&d[i]);
37 }
39 u[1]=b[1];
40 y[1]=d[1];
41 for(i=2; i<=n; i++) {
42 l[i]=a[i]/u[i−1];
43 u[i]=b[i]−l[i]*c[i−1];
44 y[i]=d[i]−l[i]*y[i−1];
45 }
46 x[n]=y[n]/u[n];
47 for(i=n−1; i>0; i−−)
48 x[i]=(y[i]−c[i]*x[i+1])/u[i];
49 cprintf(”The result:\r\n\n”);
50 for(i=n; i>=1; i−−)
51 printf(” x%d=%f\n”,i,x[i]);
52 getche();
53 printf(”\n”);
54 printf(”continue?(y/n)”);
55 play=getche();
56 } while(play==’y’||play==’Y’);
57
58 }

Jacobi 迭代法

#define MAX 100
4 #define n 3
5 #define exp 0.005
6 main() {
7 int i,j,k,m;
8 float temp,s;
9 float static a[n][n]= {{10,0,−1},{−2,10,−1},{0,−1,5}};
10 float static b[n]= {9,7,4};
11 float static x[n],B[n][n],g[n],y[n]= {0,0,0};
13 for(i=0; i<n; i++)
14 for(j=0; j<n; j++) {
15 B[i][j]=−a[i][j]/a[i][i];
16 g[i]=b[i]/a[i][i];
17 }
18 for(i=0; i<n; i++)
19 B[i][i]=0;
20 m=0;
21 do {
22 for(i=0; i<n; i++)
23 x[i]=y[i];
24 for(i=0; i<n; i++) {
25 y[i]=g[i];
26 for(j=0; j<n; j++)
27 y[i]=y[i]+B[i][j]*x[j];
28 }
29 m++;
30 printf(”\n%dth result is:”,m);
31 printf(”\nx0=%7.5f,x1=%7.5f,x2=%7.5f”,y[0],y[1],y[2]);
32 temp=0;
33 for(i=0; i<n; i++) {
34 s=fabs(y[i]−x[i]);
35 if(temp<s) temp=s;
36 }
37 printf(”\ntemp=%f”,temp);
38 } while(temp>=exp);
39 printf(”\n\nThe last result is:”);
40 for(i=0; i<n; i++)
41 printf(”\nx[%d]=%7.5f”,i,y[i]);
42 }

Gauss-Seidel 迭代法

#define MAXSIZE 50
5 void input(double a[MAXSIZE][MAXSIZE],double b[MAXSIZE],long n);
6 void output(double x[MAXSIZE],long n);
9 void main(void) {
10 double a[MAXSIZE][MAXSIZE],b[MAXSIZE],x[MAXSIZE];
11 double epsilon,e,s,oldx;
12 long n,i,j,k,maxk;
13 printf(”\n请输入原方程组的阶数：”);
14 scanf(”%ld”,&n);
15 input(a,b,n);
16 printf(”\n请输入迭代初始向量：”);
17 for(i=0; i<=n−1; i++)
18 scanf(”%lf”,&x[i]);
19 printf(”\n请输入最大迭代次数：”);
20 scanf(”%ld”,&maxk);
21 printf(”\n请输入误差上限：”);
22 scanf(”%lf”,&epsilon);
23 for(k=1; k<=maxk; k++) {
24 e=0;
25 for(i=0; i<=n−1; i++) {
26 oldx=x[i];
27 s=0;
28 for(j=0; j<=n−1; j++)
29 if(j!=i) s+=a[i][j]*x[j];
30 x[i]=(b[i]−s)/a[i][i];
31 if(e<fabs(oldx−x[i]))
32 e=fabs(oldx−x[i]);
33 }
34 if(e<epsilon) break;
35 }
36 if(k<=maxk)
37 output(x,n);
38 else
39 printf(”\n迭代次数已超过上限。”);
40 }
42 void input(double a[MAXSIZE][MAXSIZE],double b[MAXSIZE],long n) {
43 long i,j;
44 printf(”\n请输入原方程组的增广矩阵：\n”);
45 for(i=0; i<=n−1; i++) {
46 for(j=0; j<=n−1; j++)
47 scanf(”%lf”,&a[i][j]);
48 scanf(”%lf”,&b[i]);
49 }
50 }
52 void output(double x[MAXSIZE],long n) {
53 long i;
54 printf(”\n原方程组的解向量为：\n”);
55 for(i=0; i<=n−1; i++)
56 printf(” %lf\n”,x[i]);
57 }