Fortran Coder

查看: 135|回复: 4

[求助] Doolittle分解的算法实现问题

[复制链接]

13

帖子

4

主题

0

精华

入门

F 币
60 元
贡献
30 点
发表于 2024-3-30 11:41:17 | 显示全部楼层 |阅读模式
各位老师好,我在写用“Doolittle分解法解线性方程组”这一数值算法的实现程序时遇到了问题,程序给出了错误的结果,想请教各位老师这个程序在算法上有什么错误,谢谢!
(附件中是实现doolittle分解法的子程序,这是我写的交互程序的一部分,A,X,B分别代表系数矩阵,解,常数矩阵;ndim代表方阵A的维度)
Doolittle_Decomposition.f90 (2.01 KB, 下载次数: 2)

5

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
16 点
发表于 2024-3-30 11:54:16 | 显示全部楼层
你在换列的时候是不是写错写成换行了

5

帖子

1

主题

0

精华

入门

F 币
50 元
贡献
16 点
发表于 2024-3-30 12:00:10 | 显示全部楼层
[Fortran] 纯文本查看 复制代码
      module matixsol
      implicit none
      contains
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c matrix solving program
c 2024.1.18
c  type|variable | description
c   i  |mata| linear coefficient matrix
c   i  |matb| constant array
c   o  |matx| element array
c the matrixs must be full rank
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine matsolve(mata,matb,matx)
      implicit none
c input/output
      real*8 :: mata(:,:)
      real*8 :: matb(:)
      real*8 :: matx(:)
c temporary variates
      real*8,allocatable :: matl(:,:)
      real*8,allocatable :: matu(:,:)
      real*8,allocatable :: maty(:)
      integer*4 :: numb
      integer*4 :: cycl

ccccccccccccccccccccccccccccccccccc
      numb = size( matb )
      allocate( matl(numb,numb) )
      allocate( matu(numb,numb) )
      allocate( maty(numb) )
c initialize
      matl = 0.0
      matu = 0.0
      maty = 0.0

c LU decomposition of matrix
      matu(1 , 1:numb) = mata(1 , 1:numb)
      matl(1:numb , 1) = mata(1:numb , 1) / matu(1 , 1)
      do cycl = 2 , numb
        matu(cycl , cycl:numb) = mata(cycl , cycl:numb)
     &                         - matmul(matl(cycl , 1:cycl-1)
     &                         , matu(1:cycl-1 , cycl:numb))
        if (cycl .le. numb) then
          matl(cycl:numb , cycl) = (mata(cycl : numb , cycl)
     &                           - matmul(matl(cycl:numb
     &                           , 1:cycl-1) , matu( 1:cycl-1
     &                           , cycl))) / matu(cycl , cycl)
        end if
      end do
c solve Uy=b
      maty(1) = matb(1)
      do cycl = 2 , numb
        maty(cycl) = matb(cycl) - dot_product(matl(cycl , 1:cycl-1)
     &             , maty(1:cycl-1))
      end do

c solve Ux=y
      matx(numb) = maty(numb) / matu(numb , numb)               
      do cycl = numb-1 , 1 , -1
        matx(cycl) = (maty(cycl) - dot_product( matu(cycl , cycl + 1
     &             : numb),matx(cycl + 1 : numb))) / matu(cycl,cycl)
      end do
      return
      end subroutine matsolve
      end module matixsol

这是我之前写的LU的代码,里面也没有换行的操作

评分

参与人数 1F 币 +20 收起 理由
Yjc + 20 谢谢老师,我已经解决了

查看全部评分

13

帖子

4

主题

0

精华

入门

F 币
60 元
贡献
30 点
 楼主| 发表于 2024-3-30 12:19:21 | 显示全部楼层
感谢老师解答,我的想法是在LU分解中结合列选主元,把每列中最大的元素放到主元所在的位置,以增加数值稳定性

13

帖子

4

主题

0

精华

入门

F 币
60 元
贡献
30 点
 楼主| 发表于 2024-3-30 13:55:49 | 显示全部楼层
我发现问题所在了,在选列主元的时候,应该同时对矩阵A和矩阵L进行行变换,使得主元位置元素的绝对值最大。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-4-22 18:17

Powered by Tencent X3.4

© 2013-2024 Tencent

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