用C语言写一个高斯消元法解方程组的程序

我们以方程组2x1 + 6x2 - x3 = -12

5x1 - x2 +2x3 = 29

-3x1 - 4x2 + x3 = 5

为例 来说明楼主自己把方程组化为矩阵形式。以下为源代码 。

#include <stdio.h>

#include <stdlib.h>

#include <malloc.h>

#include <math.h>

int GS(int,double**,double *,double)

double **TwoArrayAlloc(int,int)

void TwoArrayFree(double **)

int main(void)

{

int i,n

double ep,**a,*b

n = 3

ep = 1e-4

a = TwoArrayAlloc(n,n)

b = (double *)calloc(n,sizeof(double))

if(b == NULL)

{

printf("memory get error\n")

exit(1)

}

a[0][0]= 2a[0][1]= 6a[0][2]= -1

a[1][0]= 5a[1][1]=-1a[1][2]= 2

a[2][0]=-3a[2][1]=-4a[2][2]= 1

b[0] = -12b[1] = 29 b[2] = 5

if(!GS(n,a,b,ep))

{

printf("can't solve it with GS elimination\n")

exit(0)

}

printf("The solution of equations is as follows：\n")

for(i=0i<3i++)

{

printf("x%d = %.2f\n",i,b[i])

}

TwoArrayFree(a)

free(b)

return 0

}

int GS(n,a,b,ep)

int n

double **a

double *b

double ep

{

int i,j,k,l

double t

for(k=1k<=nk++)

{

for(l=kl<=nl++)

if(fabs(a[l-1][k-1])>ep)

break

else if(l==n)

return(0)

if(l!=k)

{

for(j=kj<=nj++)

{

t = a[k-1][j-1]

a[k-1][j-1] =a[l-1][j-1]

a[l-1][j-1] =t

}

t=b[k-1]

b[k-1]=b[l-1]

b[l-1]=t

}

t=1/a[k-1][k-1]

for(j=k+1j<=nj++)

a[k-1][j-1]=t*a[k-1][j-1]

b[k-1]*=t

for(i=k+1i<=ni++)

{

for(j=k+1j<=nj++)

a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1]

b[i-1]-=a[i-1][k-1]*b[k-1]

}

}

for(i=n-1i>=1i--)

for(j=i+1j<=nj++)

b[i-1]-=a[i-1][j-1]*b[j-1]

return(1)

}

double **TwoArrayAlloc(int r,int c)

{

double *x,**y

int n

x=(double *)calloc(r*c,sizeof(double))

y=(double **)calloc(r,sizeof(double*))

for(n=0n<=r-1++n)

{

y[n]=&x[c*n]

}

return y

}

void TwoArrayFree(double **x)

{

free(x[0])

free(x)

}

#include <stdio.h>

#include <iostream.h>

#include <math.h>

#include<iomanip.h>

class Angle

{

public:

double degree,cent,second,Hudu,seconds

//

构造函数

Angle(double _degree,double _cent,double _second)

{

degree=_degree

cent=_cent

second=_second

seconds=(degree*3600+cent*60+second)

Hudu=(degree*3600+cent*60+second)*2*3.1415926535/(360*60*60)

}

Angle()

{

}

SToD()

{

second=seconds-int(seconds/60)*60

degree=int((seconds-second)/3600)

cent=int((seconds-second-degree*3600)/60)

}

~Angle()

{

}

}

void main()

{

//

正算过程

//

定义变量

Angle B(51,38,43.9023),L(126,2,13.1360),_B,_L

double L0,l,N,a0,a4,a6,a3,a5,x,y,rou

rou=(360*60*60)/(2*3.1415926535)

if(int(L.degree)%6!=0)

{

L0=6*(int(L.degree)/6+1)-3+6

}

else

L0=6*(int(L.degree)/6)-3+6

l=L.Hudu-L0*3600*2*3.1415926535/(360*60*60)

N=6399698.902-(21562.267-(108.973-0.612*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(

B.Hudu))*cos(B.Hudu)*cos(B.Hudu)

a0=32140.404-(135.3302-(0.7092-0.0040*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.

Hudu))*cos(B.Hudu)*cos(B.Hudu)

a4=(0.25+0.00252*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hudu)-0.04166

a6=(0.166*cos(B.Hudu)*cos(B.Hudu)-0.084)*cos(B.Hudu)*cos(B.Hudu)

a3=(0.3333333+0.001123*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hudu)-0.166666

7

a5=0.0083-(0.1667-(0.1968+0.0040*cos(B.Hudu)*cos(B.Hudu))*cos(B.Hudu)*cos(B.Hudu)

)*cos(B.Hudu)*cos(B.Hudu)

x=6367558.4969*B.Hudu-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B.Hudu)*cos(B.Hudu)

y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B.Hudu)

//结果输出

cout<<"x="<<x<<endl

cout<<"y="<<y<<endl

// 反算过程

double b,Bf,Nf,Z,b2,b3,b4,b5,l1

b=(x/6367558.4969)

Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b))*cos(b)*cos(b))*0.000

0000001*sin(b)*cos(b)

Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf)*c

os(Bf)

Z=y/(Nf*cos(Bf))

b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf)

b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf)

b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf)

b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf)

_B.seconds=Bf*rou-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2*rou

l1=(1-(b3-b5*Z*Z)*Z*Z)*Z*rou

_L.seconds=L0*3600+l1

_B.SToD()

_L.SToD()

---