一物理系统可用下列线性方程组来表示：

一物理系统可用下列线性方程组来表示：
从文件中读入m1、m2和θ的值,求a1、a2、N1 和N2的值.其中g取9.8,输入θ时以角度为单位.
要求：
（1）分别用两种方法（例如高斯消去法、矩阵求逆法、三角分解法、追赶法等）,定义求解线性方程组Ax=b的子程序,要求该子程序能求解任意线性方程组.
（2）在主程序中分别调用上面定义的两个子程序,并对求解结果进行对比分析.
（3）绘制以上两个方法所求得的方程解的数据分布图.
谢谢(立模块3个,子程序,文件调用,高斯消去)用FORTRAN语言

SUBROUTINE GAUSS(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) !将对角线上元素变为1
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 !化为每行主对角线左边的数为0的矩阵
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 MAIN1
DIMENSION A(4,4),B(4),X(4),JS(4)
DOUBLE PRECISION A,B,X
REAL M1,M2,C
OPEN(1,FILE="zhouxiao.TXT")
READ(1,*)M1,M2,C
CLOSE(1)
N=4
PRINT*,M1,M2,C
A(1,1)=M1*COS(3.14159*C/180)
A(1,2)=-M1
A(1,3)=-SIN(3.14159*C/180)
A(1,4)=0
A(2,1)=M1*SIN(3.14159*C/180)
A(2,2)=0
A(2,3)=COS(3.14159*C/180)
A(2,4)=0
A(3,1)=0
A(3,2)=M2
A(3,3)=-SIN(3.14159*C/180)
A(3,4)=0
A(4,1)=0
A(4,2)=0
A(4,3)=-COS(3.14159*C/180)
A(4,4)=1
B(1)=0
B(2)=M1*9.8
B(3)=0
B(4)=M2*9.8
CALL GAUSS(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