用直接消去法解方程组的程序如何编写？(Fortran程序）

!高纤册斯消去法

subroutine agaus(a,b,n,x,l,js)

dimension a(n,n),x(n),b(n),js(n)

double precision a,b,x,t

l=1 !逻辑变量

do k=1,n-1

d=0.0

do i=k,n

do j=k,n

if (abs(a(i,j))>d) then

d=abs(a(i,j))

js(k)=j

is=i

end if

end do

end do !把行绝对值最大的元素换到粗蠢主元位毁凳宏置

if (d+1.0==1.0) then

l=0

else !最大元素为0无解

if(js(k)/=k) then

do i=1,n

t=a(i,k)

a(i,k)=a(i,js(k))

a(i,js(k))=t

end do !最大元素不在K行，K行

end if

if(is/=k) then

do j=k,n

t=a(k,j)

a(k,j)=a(is,j)

a(is,j)=t !交换到K列

end do

t=b(k)

b(k)=b(is)

b(is)=t

end if !最大元素在主对角线上

end if !消去

if (l==0) then

write(*,100)

return

end if

do j=k+1,n

a(k,j)=a(k,j)/a(k,k)

end do

b(k)=b(k)/a(k,k) !求三角矩阵

do i=k+1,n

do j=k+1,n

a(i,j)=a(i,j)-a(i,k)*a(k,j)

end do

b(i)=b(i)-a(i,k)*b(k)

end do

end do

if (abs(a(n,n))+1.0==1.0) then

l=0

write(*,100)

return

end if

x(n)=b(n)/a(n,n)

do i=n-1,1,-1

t=0.0

do j=i+1,n

t=t+a(i,j)*x(j)

end do

x(i)=b(i)-t

end do

100 format(1x,'fail')

js(n)=n

do k=n,1,-1

if (js(k)/=k) then

t=x(k)

x(k)=x(js(k))

x(js(k))=t

end if

end do

return

end

program main

dimension a(4,4),b(4),x(4),js(4)

double precision a,b,x

real m1,m2,j

open(1,file="laiyi.txt")

read(1,*)m1,m2,j

close(1)

n=4

print*,m1,m2,j

a(1,1)=m1*cos(3.14159*j/180)

a(1,2)=-m1

a(1,3)=-sin(3.14159*j/180)

a(1,4)=0

a(2,1)=m1*sin(3.14159*j/180)

a(2,2)=0

a(2,3)=cos(3.14159*j/180)

a(2,4)=0

a(3,1)=0

a(3,2)=m2

a(3,3)=-sin(3.14159*j/180)

a(3,4)=0

a(4,1)=0

a(4,2)=0

a(4,3)=-cos(3.14159*j/180)

a(4,4)=1

b(1)=0

b(2)=m1*9.8

b(3)=0

b(4)=m2*9.8

call agaus(a,b,n,x,l,js)

if (l/=0) then

write(*,*)"a1=",x(1),"a2=",x(2) ,"n1=",x(3),"n2=",x(4)

end if

end

!逆矩阵求解

SUBROUTINE qiuni(A,N,L,IS,JS)

DIMENSION A(N,N),IS(N),JS(N)

DOUBLE PRECISION A,T,D

L=1

DO K=1,N

D=0.0

DO I=K,N

DO J=K,N

IF(ABS(A(I,J)).GT.D) THEN !把最大的元素给D

D=ABS(A(I,J))

IS(K)=I

JS(K)=J

END IF

END DO

END DO

IF (D+1.0.EQ.1.0)THEN

L=0

WRITE(*,200)

RETURN

END IF

200 FORMAT(1X,'ERR**NOT INV')

DO J=1,N

T=A(K,J)

A(K,J)=A(IS(K),J)

A(IS(K),J)=T

END DO

DO I=1,N

T=A(I,K)

A(I,K)=A(I,JS(K))

A(I,JS(K))=T

END DO

A(K,K)=1/A(K,K)

DO J=1,N

IF(J.NE.K)THEN

A(K,J)=A(K,J)*A(K,K)

END IF

END DO

DO I=1,N

IF(I.NE.K)THEN

DO J=1,N

IF(J.NE.K)THEN

A(I,J)=A(I,J)-A(I,K)*A(K,J)

END IF

END DO

END IF

END DO

DO I=1,N

IF(I.NE.K)THEN

A(I,K)=-A(I,K)*A(K,K)

END IF

END DO

END DO

DO K=N,1,-1

DO J=1,N

T=A(K,J)

A(K,J)=A(JS(K),J)

A(JS(K),J)=T

END DO

DO I=1,N

T=A(I,K)

A(I,K)=A(I,IS(K))

A(I,IS(K))=T

END DO

END DO

RETURN

END

SUBROUTINE BRMUL(A,B,N,C)

DIMENSION A(N,N),B(N),C(N)

DOUBLE PRECISION A,B,C

DO I=1,N

DO J=1,N

C(I)=0.0

DO L=1,N

C(I)=C(I)+A(I,L)*B(L)

END DO

END DO

END DO

RETURN

END

program main

DIMENSION A(4,4),B(4,1),C(4,1),IS(4),JS(4)

DOUBLE PRECISION A,B,C

REAL M1,M2,JD

OPEN(1,FILE='LAIYI.TXT')

READ(1,*) M1,M2,JD

PRINT*,M1,M2,JD

CLOSE(1)

A(1,1)=M1*COS(3.14*JD/180)

A(1,2)=-M1

A(1,3)=-SIN(3.14*JD/180)

A(1,4)=0

A(2,1)=M1*SIN(3.14*JD/180)

A(2,2)=0

A(2,3)=COS(3.14*JD/180)

A(2,4)=0

A(3,1)=0

A(3,2)=M2

A(3,3)=-SIN(3.14*JD/180)

A(3,4)=0

A(4,1)=0

A(4,2)=0

A(4,3)=-COS(3.14*JD/180)

A(4,4)=1

B(1,1)=0

B(2,1)=M1*9.8

B(3,1)=0

B(4,1)=M2*9.8

CALL QIUNI(A,4,L,IS,JS)

CALL BRMUL(A,B,4,C)

WRITE(*,*) (C(I,1),I=1,4)

END

画图

USE MSFLIB

INTEGER status

TYPE(xycoord) xy

status=SETCOLORRGB(#FFFFFF)

status1=SETCOLORRGB(#0000FF)

OPEN(1,FILE="G.TXT")

READ(1,*) G1,G2,G3,G4

OPEN(2,FILE="N.TXT")

READ(2,*) N1,N2,N3,N4

CALL MOVETO(INT(20),INT(20),xy)

status=LINETO(INT(40),INT(G1))

status=LINETO(INT(80),INT(G2))

status=LINETO(INT(120),INT(G3))

status=LINETO(INT(160),INT(G4))

CALL SETLINESTYLE(#FF00)

CALL MOVETO(INT(20),INT(20),xy)

status1=LINETO(INT(40),INT(N1))

status1=LINETO(INT(80),INT(N2))

status1=LINETO(INT(120),INT(N3))

status1=LINETO(INT(160),INT(N4))

READ(*,*)

END

如果对您有帮助，请记得采纳为满意答案，谢谢！祝您生活愉快！

#include<iostream>

#include<cmath>

using namespace std

#define MAX 50

void input(double a[MAX][MAX+1],int n)

{

cout<<"输入原方程组的增广矩阵"<<endl

for(int i=0i<ni++)

for(int j=0j<n+1j++)

cin>>a[i][j]

}

void output(double x[],int n)

{

cout<<"Gauss 消去法得到的原方程组的解为"<<endl

for(int k=0k<nk++)

cout<<x[k]<<" "

}

int main()

{

double a[MAX][MAX+1],x[MAX],sum,max,t

int n,i,j,k,max_i

cout<逗仔<"输入原方程组的阶"<<endl cin>>n

input(a,n)

for(k=0k<n-1k++)//选主元素

{ max=a[k][k]

max_i=k

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

if(fabs(a[i][k])>fabs(max))

{

max=a[i][k]

max_i=i

}

if(max==0)

break

if(max_i!=k)//交换两行脊迹

for(j=kj<n+1j++)

{

t=a[k][j]

a[k][j]=a[max_i][j]

a[max_i][j]=t

}

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

{

a[i][k]=a[i][k]/-a[k][k]

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

a[i][j]=a[i][j]+a[i][k]*a[k][j]

}//消元

}

if(max==0)cout<<"原方程组无解"<<endl

else

{

for(k=n-1k>=0k--)

{

sum=0

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

sum=sum+a[k][j]*x[j]

x[k]=(a[k][n]-sum)/a[k][k]

}//回代

output(x,n)

cout<樱指并<endl

}

return 0

}

A=[6,-2,2,412,-8,4,103,-13,3,3-6,4,2,18]

b=[0-10-39-16]

B=[A b]

ra=rank(A)

rb=rank(B)

n=length(b)

X=zeros(n,1)C=zeros(1,n+1)

if ra>rb

disp('B的秩大于A的秩，方程组无解！')

return

elseif ra==rb &rb==n

disp('B和A的秩均等型握于n,方程组有唯一解，用高斯列主元消去法求解过程如下：')

for p=1:n-1

disp('卜困庆p=')

disp(p)

[Y,j]=max(abs(B(p:n,p)))

disp('j=')

disp(j)

C=B(p,:)

B(p,:)= B(j+p-1,:)B(j+p-1,:)=C

disp(B)

for k=p+1:n

m= B(k,p)/ B(p,p)

B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1)

end

disp(B)

end

b=B(1:n,n+1)A=B(1:n,1:n)X(n)=b(n)/A(n,n)

for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q)

end

disp(X)

else

disp('B和A的秩相等切尺梁小于n,方程组有无穷解！')

end

---