Fortran Coder

查看: 241|回复: 7

[线性代数] 如何用Fortran求解5对角矩阵?

[复制链接]

4

帖子

2

主题

0

精华

入门

F 币
45 元
贡献
22 点
发表于 2020-7-2 08:56:48 | 显示全部楼层 |阅读模式
如何用Fortran求解5对角矩阵啊,有没有大佬有程序贴出来让我学习一下,谢谢!
回复

使用道具 举报

1538

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1417 元
贡献
936 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2020-7-2 10:57:40 | 显示全部楼层
矩阵不存在“求解”问题。它本身只是数据表达的一种方式,不是一个“待解的问题”。
你需要描述你的具体问题,比如:
求解系数矩阵为5阶对角矩阵的线性问题 Ax=b
或者,求解5阶对角矩阵的逆矩阵,或者某种分解

4

帖子

2

主题

0

精华

入门

F 币
45 元
贡献
22 点
 楼主| 发表于 2020-7-2 15:02:40 | 显示全部楼层
fcode 发表于 2020-7-2 10:57
矩阵不存在“求解”问题。它本身只是数据表达的一种方式,不是一个“待解的问题”。
你需要描述你的具体问 ...

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


546

帖子

1

主题

0

精华

大宗师

F 币
2702 元
贡献
1721 点
发表于 2020-7-2 22:26:42 | 显示全部楼层
你的需求是解方程,不管是三对角还是五对角,甚至任意矩阵,都可以用  通用函数 来求解,找一个解方程的函数就是,无外乎效率的差异。

116

帖子

8

主题

0

精华

版主

F 币
557 元
贡献
336 点

贡献勋章管理勋章帅哥勋章元老勋章星光勋章规矩勋章

发表于 2020-7-5 00:19:17 | 显示全部楼层

RE: 如何用Fortran求解5对角矩阵?

[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                                        

請參考:

tri-penta.pdf

1.88 MB, 下载次数: 2

4

帖子

2

主题

0

精华

入门

F 币
45 元
贡献
22 点
 楼主| 发表于 2020-7-5 08:39:33 | 显示全部楼层
chiangtp 发表于 2020-7-5 00:19
[mw_shl_code=fortran,false]!---------------------------------------------------------------------
! ...

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

74

帖子

18

主题

1

精华

大师

F 币
743 元
贡献
446 点
发表于 2020-8-6 21:54:44 | 显示全部楼层
本帖最后由 weixing1531 于 2020-8-6 22:00 编辑
Pluto 发表于 2020-7-2 15:02
如图所示。我问的其实是数值传热学中的一个问题,书中给出了求解三对角阵算法的过程,现在我要求解五对角 ...

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

160

帖子

14

主题

0

精华

大宗师

F 币
3921 元
贡献
774 点
发表于 5 天前 | 显示全部楼层
为什么这种事情还要去研究它的细节……
好多现成的库函数来解决这种问题吧……
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-8-14 22:07

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表