如何用Fortran求解5对角矩阵？

Pluto

如何用Fortran求解5对角矩阵啊，有没有大佬有程序贴出来让我学习一下，谢谢！

fcode

矩阵不存在“求解”问题。它本身只是数据表达的一种方式，不是一个“待解的问题”。
你需要描述你的具体问题，比如：
求解系数矩阵为5阶对角矩阵的线性问题 Ax=b
或者，求解5阶对角矩阵的逆矩阵，或者某种分解？

Pluto

fcode 发表于 2020-7-2 10:57
矩阵不存在“求解”问题。它本身只是数据表达的一种方式，不是一个“待解的问题”。
你需要描述你的具体问 ...

如图所示。我问的其实是数值传热学中的一个问题，书中给出了求解三对角阵算法的过程，现在我要求解五对角阵算法的过程，也就是求
,相当于求二维的情况。这种应该怎么算

li913

你的需求是解方程，不管是三对角还是五对角，甚至任意矩阵，都可以用  通用函数 来求解，找一个解方程的函数就是，无外乎效率的差异。

chiangtp

[Fortran] 纯文本查看 复制代码

查看代码

打印?

!---------------------------------------------------------------------

! Routine to solve a diagonally dominant scalar tridiagonal system

! by compact Crout LU method.

!

!    A(I)*x(I-1) + B(I)*x(I) + C(I)*x(I+1) = D(I)

!       for 1 <= I <= N, with A(1) = C(N) = 0

!

!   A   : the subdiagonal coefficients, size N,

!           A(1) is not used.

!   B   : the diagonal coefficients, size N.

!   C   : the above-diagonal coefficients, size N,

!           C(N) is not used.

!   D   : the right-hand-side coefficients, size N,

!           will contain the solution x on exit

!   Ier : Error flag,

!           =  0, normal return

!           = -1, zero-sized of A, no problem to solve

!           = -2, array B, C, or, D not conformable with A

!---------------------------------------------------------------------

SUBROUTINE Axb_Tri_Scalar( A, B, C, D, Ier )

IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN   ) :: A

REAL, DIMENSION(:), INTENT(INOUT) :: B, C, D

INTEGER,            INTENT(  OUT) :: Ier

INTEGER :: N, I

!-------------------------------------

N = SIZE( A )

IF( N == 0 ) THEN

Ier = -1

RETURN

END IF

IF( ( N /= SIZE(B) ) .OR. &

( N /= SIZE(C) ) .OR. &

( N /= SIZE(D) ) ) THEN

Ier = -2

RETURN

END IF

!-----------------

C(N) = 0.0

C(1) = C(1) / B(1)

D(1) = D(1) / B(1)

DO I = 2, N

B(I) =   B(I) - A(I)*C(I-1)

D(I) = ( D(I) - A(I)*D(I-1) ) / B(I)

C(I) = C(I) / B(I)

END DO

! Backward Substitution

DO I = (N-1), 1, -1

D(I) = D(I) - C(I)*D(I+1)

END DO

Ier = 0

END SUBROUTINE Axb_Tri_Scalar

!---------------------------------------------------------------------

! Routine to solve a diagonally dominant scalar pentadiagonal system

! by compact Crout LU method.

!

! A(I)*x(I-2)+B(I)*x(I-1)+C(I)*x(I)+D(I)*x(I+1)+E(I)*x(I+2)=F(I)

!    for 1 <= I <= N, with A(1)=A(2)=B(1)=D(N)=E(N-1)=E(N)=0

!

!   A   : the 1st subdiagonal coefficients, size N,

!           A(1) and A(2) are not used.

!   B   : the 2nd subdiagonal coefficients, size N,

!           B(1) is not used.

!   C   : the diagonal coefficients, size N.

!   D   : the 1st above-diagonal coefficients, size N,

!           D(N) is not used.

!   E   : the 2nd above-diagonal coefficients, size N,

!           E(N) and E(N-1) are not used.

!   F   : the right-hand-side coefficients, size N,

!           will contain the solution x on exit

!   Ier : Error flag,

!           =  0, normal return

!           = -1, zero-sized of A, no problem to solve

!           = -2, array B,C,D,E, or, F not conformable with A

!---------------------------------------------------------------------

SUBROUTINE Axb_Penta_Scalar( A, B, C, D, E, F, Ier )

IMPLICIT NONE

REAL, DIMENSION(:), INTENT(IN   ) :: A

REAL, DIMENSION(:), INTENT(INOUT) :: B, C, D, E, F

INTEGER,            INTENT(  OUT) :: Ier

INTEGER :: N, N1, N2, N3, I, I1, I2

!-------------------------------------

N = SIZE(A)

IF( N == 0 ) THEN

Ier = -1

RETURN

END IF

IF( ( N /= SIZE(B) ) .OR. &

( N /= SIZE(C) ) .OR. &

( N /= SIZE(D) ) .OR. &

( N /= SIZE(E) ) .OR. &

( N /= SIZE(F) ) ) THEN

Ier = -2

RETURN

END IF

!-----------------

D(1) = D(1) / C(1)

E(1) = E(1) / C(1)

F(1) = F(1) / C(1)

C(2) =   C(2) - B(2)*D(1)

D(2) = ( D(2) - B(2)*E(1) ) / C(2)

E(2) =   E(2)               / C(2)

F(2) = ( F(2) - B(2)*F(1) ) / C(2)

DO I = 3, (N-2)

I1 = I-1

I2 = I-2

B(I) =   B(I)              - A(I)*D(I2)

C(I) =   C(I) - B(I)*D(I1) - A(I)*E(I2)

D(I) = ( D(I) - B(I)*E(I1)              ) / C(I)

E(I) =   E(I)                             / C(I)

F(I) = ( F(I) - B(I)*F(I1) - A(I)*F(I2) ) / C(I)

END DO

N1 = N-1

N2 = N-2

N3 = N-3

B(N1) =   B(N1)               - A(N1)*D(N3)

C(N1) =   C(N1) - B(N1)*D(N2) - A(N1)*E(N3)

D(N1) = ( D(N1) - B(N1)*E(N2)               ) / C(N1)

F(N1) = ( F(N1) - B(N1)*F(N2) - A(N1)*F(N3) ) / C(N1)

N1 = N-1

N2 = N-2

B(N) =   B(N)              - A(N)*D(N2)

C(N) =   C(N) - B(N)*D(N1) - A(N)*E(N2)

F(N) = ( F(N) - B(N)*F(N1) - A(N)*F(N2) ) /C(N)

!----------------------

! Backward Substitution

!----------------------

I = N-1

F(I) = F(I) - D(I)*F(I+1)

DO I = (N-2), 1, -1

F(I) = F(I) - D(I)*F(I+1) - E(I)*F(I+2)

END DO

Ier = 0

END SUBROUTINE Axb_Penta_Scalar

請參考:

Pluto

chiangtp 发表于 2020-7-5 00:19
[mw_shl_code=fortran,false]!---------------------------------------------------------------------
! ...

谢谢，我去研究研究这个程序

weixing1531

Pluto 发表于 2020-7-2 15:02
如图所示。我问的其实是数值传热学中的一个问题，书中给出了求解三对角阵算法的过程，现在我要求解五对角 ...

图片上这本书是陶文铨院士所著的《数值传热学》
西安交大网站上有SIMPLE算法的源代码http://nht.xjtu.edu.cn/xnfzsyjxzx/zyxz.htm
思想：五对角转换成三对角   先TDMA后ADI（交替方向隐格式）

liudy02

为什么这种事情还要去研究它的细节……
好多现成的库函数来解决这种问题吧……