!高斯消去法
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=00
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+10==10) 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))+10==10) then
l=0
write(,100)
return
end if
x(n)=b(n)/a(n,n)
do i=n-1,1,-1
t=00
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="laiyitxt")
read(1,)m1,m2,j
close(1)
n=4
print,m1,m2,j
a(1,1)=m1cos(314159j/180)
a(1,2)=-m1
a(1,3)=-sin(314159j/180)
a(1,4)=0
a(2,1)=m1sin(314159j/180)
a(2,2)=0
a(2,3)=cos(314159j/180)
a(2,4)=0
a(3,1)=0
a(3,2)=m2
a(3,3)=-sin(314159j/180)
a(3,4)=0
a(4,1)=0
a(4,2)=0
a(4,3)=-cos(314159j/180)
a(4,4)=1
b(1)=0
b(2)=m198
b(3)=0
b(4)=m298
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=00
DO I=K,N
DO J=K,N
IF(ABS(A(I,J))GTD) THEN !把最大的元素给D
D=ABS(A(I,J))
IS(K)=I
JS(K)=J
END IF
END DO
END DO
IF (D+10EQ10)THEN
L=0
WRITE(,200)
RETURN
END IF
200 FORMAT(1X,'ERRNOT 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(JNEK)THEN
A(K,J)=A(K,J)A(K,K)
END IF
END DO
DO I=1,N
IF(INEK)THEN
DO J=1,N
IF(JNEK)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(INEK)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)=00
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='LAIYITXT')
READ(1,) M1,M2,JD
PRINT,M1,M2,JD
CLOSE(1)
A(1,1)=M1COS(314JD/180)
A(1,2)=-M1
A(1,3)=-SIN(314JD/180)
A(1,4)=0
A(2,1)=M1SIN(314JD/180)
A(2,2)=0
A(2,3)=COS(314JD/180)
A(2,4)=0
A(3,1)=0
A(3,2)=M2
A(3,3)=-SIN(314JD/180)
A(3,4)=0
A(4,1)=0
A(4,2)=0
A(4,3)=-COS(314JD/180)
A(4,4)=1
B(1,1)=0
B(2,1)=M198
B(3,1)=0
B(4,1)=M298
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="GTXT")
READ(1,) G1,G2,G3,G4
OPEN(2,FILE="NTXT")
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
点除与矩阵除法:
在书写程序的时候,点乘和矩阵乘法写错的时候再进行程序调适的
时候MATLAB会返回错误说明。
但是对于点除容易出现问题,下面以一个简单的例子说明这个问题:
比如我们要计算:
A = [1,1];
B = [2,1];
C = A/B;
上面的程序我们计算的是A与B的点除。但是由于疏忽而把点除“/”
以上就是关于fortran:用高斯消去法,矩阵求逆法、三角分解法、追赶法中的两种解线性方程组全部的内容,包括:fortran:用高斯消去法,矩阵求逆法、三角分解法、追赶法中的两种解线性方程组、为什么用追赶法和矩阵除法算的结果不一样、等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)