一物理系统可用下列线性方程组来表示：从文件中读入m1、m2和θ的值,求a1、a2、N1和N2的值.其中g取9.8,输入θ时以角度为单位.要求：（1）分别用两种方法（例如高斯消去法、矩阵求逆法、三
一物理系统可用下列线性方程组来表示：
从文件中读入m1、m2和θ的值,求a1、a2、N1和N2的值.其中g取9.8,输入θ时以角度为单位.
要求：
（1）分别用两种方法（例如高斯消去法、矩阵求逆法、三角分解法、追赶法等）,定义求解线性方程组Ax=b的子程序,要求该子程序能求解任意线性方程组.
（2）在主程序中分别调用上面定义的两个子程序,并对求解结果进行对比分析.
（3）绘制以上两个方法所求得的方程解的数据分布图.
谢谢(立模块3个,子程序,文件调用,高斯消去)用FORTRAN语言
SUBROUTINEGAUSS(A,B,N,X,L,JS)
DIMENSIONA(N,N),X(N),B(N),JS(N)
DOUBLEPRECISIONA,B,X,T
L=1
DOK=1,N-1
D=0.0
DOI=K,N
DOJ=K,N
IF(ABS(A(I,J))>D)THEN
D=ABS(A(I,J))
JS(K)=J
IS=I
ENDIF
ENDDO
ENDDO!把行绝对值最大的元素换到主元位置
IF(D+1.0==1.0)THEN
L=0
ELSE!最大元素为0无解
IF(JS(K)/=K)THEN
DOI=1,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
ENDDO!最大元素不在K行交换到K行
ENDIF
IF(IS/=K)THEN
DOJ=K,N
T=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=T!将列最大元素交换到K列
ENDDO
T=B(K)
B(K)=B(IS)
B(IS)=T
ENDIF!最大元素在主对角线上
ENDIF
IF(L==0)THEN
WRITE(*,100)
RETURN
ENDIF
DOJ=K+1,N
A(K,J)=A(K,J)/A(K,K)!将对角线上元素变为1
ENDDO
B(K)=B(K)/A(K,K)!求三角矩阵
DOI=K+1,N
DOJ=K+1,N
A(I,J)=A(I,J)-A(I,K)*A(K,J)
ENDDO!化为每行主对角线左边的数为0的矩阵
B(I)=B(I)-A(I,K)*B(K)
ENDDO
ENDDO
IF(ABS(A(N,N))+1.0==1.0)THEN
L=0
WRITE(*,100)
RETURN
ENDIF
X(N)=B(N)/A(N,N)
DOI=N-1,1,-1
T=0.0
DOJ=I+1,N
T=T+A(I,J)*X(J)
ENDDO
X(I)=B(I)-T
ENDDO
100FORMAT(1X,'FAIL')
JS(N)=N
DOK=N,1,-1
IF(JS(K)/=K)THEN!上三角阵的回代
T=X(K)
X(K)=X(JS(K))
X(JS(K))=T
ENDIF
ENDDO
RETURN
END
PROGRAMMAIN1
DIMENSIONA(4,4),B(4),X(4),JS(4)
DOUBLEPRECISIONA,B,X
REALM1,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
CALLGAUSS(A,B,N,X,L,JS)
IF(L/=0)THEN
WRITE(*,*)”高斯消去法解”"A1=",X(1),"A2=",X(2),"N1=",X(3),"N2=",X(4)
ENDIF
END