fortran:用高斯消去法，矩阵求逆法、三角分解法、追赶法中的两种解线性方程组

发布网友发布时间：2022-04-29 03:08

共5个回答

热心网友时间：2023-10-08 21:15

!高斯消去法
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

热心网友时间：2023-10-08 21:15

在MKL或是IMSL里有现成的算法，您直接调用就行了。求逆是用.i.（在IMSL中）。

热心网友时间：2023-10-08 21:16

integer,parameter :: N=4
real(8),dimension(4,4) :: a
read*,a
do k=1,N
a(k,k)=1.0/a(k,k)
do j=1,N
if(j/=k) a(j,k)=a(k,k)*a(j,k)
end do
do i=1,N
do j=1,N
if(i/=k.and.j/=k) a(j,i)=a(j,i)-a(k,i)*a(j,k)
end do
if(i/=k) a(k,i)=-a(k,i)*a(k,k)
end do
end do
do i=1,N
write(*,'(4f12.6)') a(i,:)
end do
end

热心网友时间：2023-10-08 21:16

求逆矩阵 integer,parameter :: N=4
real(8),dimension(4,4) :: a
read*,a
do k=1,N
a(k,k)=1.0/a(k,k)
do j=1,N
if(j/=k) a(j,k)=a(k,k)*a(j,k)
end do
do i=1,N
do j=1,N
if(i/=k.and.j/=k) a(j,i)=a(j,i)-a(k,i)*a(j,k)
end do
if(i/=k) a(k,i)=-a(k,i)*a(k,k)
end do
end do
do i=1,N
write(*,'(4f12.6)') a(i,:)
end do
end

热心网友时间：2023-10-08 21:17

你是中南大学土建的吗？

热心网友时间：2023-10-08 21:15

在MKL或是IMSL里有现成的算法，您直接调用就行了。求逆是用.i.（在IMSL中）。

热心网友时间：2023-10-08 21:16

热心网友时间：2023-10-08 21:17

你是中南大学土建的吗？

热心网友时间：2023-10-08 21:15

在MKL或是IMSL里有现成的算法，您直接调用就行了。求逆是用.i.（在IMSL中）。

热心网友时间：2023-10-08 21:16

热心网友时间：2023-10-08 21:17

你是中南大学土建的吗？