Question

fortran编程题 线性方程组求解问题

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

Answer 1

implicit none
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)).GT.D)THEN
D=ABS(A(I,J))
JS(K)=J
IS=I
ENDIF
END DO
IF(D+1.0.EQ.1.0)THEN
L=0
ELSE
IF(JS(K).NE.k)THEN
DO I=K,N
T=A(I,K)
A(I,K)=A(I,JS(K))
A(I,JS(K))=T
END DO
ENDIF
IF(IS.NE.K)THEN
DO J=K,N
T=A(K,J)
A(K,J)=A(IS,J)
A(IS,J)=T
END DO
T=B(K)
B(K)=B(IS)
B(IS)=T
ENDIF
ENDIF
IF(L.EQ.0)THEN
WRITE(*,100)
ENDIF
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(J,K)
END DO
b(i)=b(i)-A(i,k)*b(k)
END DO
END DO
if(abs(A(n,n))+1.0.eq.1.0)then
L=0
WRITE(*,100)
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,'FALL')
JS(N)=N
DO K=N,1,-1
IF(JS(K).NE.K)THEN
T=X(K)
X(K)=X(JS(K))
X(JS(K))=T
END If
END DO
END DO