Fortran Coder

查看: 4906|回复: 3

[线性代数] 怎么就是得不出正确结果呢?

[复制链接]
回帖奖励 10 元F 币 回复本帖可获得 10 元F 币奖励! 每人限 1 次(中奖概率 50%)

35

帖子

11

主题

0

精华

熟手

F 币
129 元
贡献
118 点
发表于 2017-2-21 14:29:56 | 显示全部楼层 |阅读模式
小弟用这个高斯主元消去法解一个方程组,代码仔仔细细的看了一遍,没发现任何错误,可怎么就是结果不正确呢,请高手指点一下,很急啊,多谢多谢!!!
[Fortran] 纯文本查看 复制代码
program S_equation
implicit none
integer::i,j,k,imax,t
real::max,n

!用矩阵(实型动态数组)将线性方程组表示出来!

real,dimension(:,:),allocatable::a,m
real,dimension(:),allocatable::x
write(*,*)'请输入线性方程组的元数n:'
read(*,*)t

!给动态数组分配内存!

allocate(a(t,t+1),m(t,t+1),x(t))
write(*,*)'请依次按行输入线性方程组的增广矩阵a:'
read(*,*)((a(i,j),j=1,t+1),i=1,t)

    do k=1,t-1

    !选取每列的最大列元素!

    max=abs(a(k,k))
    imax=k
    do i=k+1,t
        if(abs(a(i,k)).gt.max)then
        max=abs(a(i,k))
        imax=i
        end if
    end do

    !将最大列元素所在的行与第K行进行交换!

    do j=k,t+1
        m(k,j)=a(k,j)
        a(k,j)=a(imax,j)
        a(imax,j)=m(k,i)
    end do

    !对方程组按X1,X2,...,Xt的顺序进行依次消元,将增广矩阵化为上三角矩阵!

    do i=k+1,t
        m(i,k)=a(i,k)/a(k,k)
        do j=k+1,t+1
            a(i,j)=a(i,j)-m(i,k)*a(k,j)
        end do
    end do
    end do

      !先算出Xt!

    x(t)=a(t,t+1)/a(t,t)

    !回带!

    do k=t-1,1,-1
        n=0
        do j=t,k+1,-1
            n=n+a(k,j)*x(j)
        end do
        x(k)=(a(k,t+1)-n)/a(k,k)
    end do

    write(*,*)'线性方程组的解为:'
    do i=1,t
        write(*,*)x(i)
    end do
    read(*,*)
    end program S_equation


1487658269(1).jpg

1958

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1339 元
贡献
565 点

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

发表于 2017-2-21 22:05:09 | 显示全部楼层
a(imax,j)=m(k,i)
改为
a(imax,j)=m(k,j)

你的代码可以进行很大的优化:
[Fortran] 纯文本查看 复制代码
program S_equation  implicit none
  integer::i,k,imax,t
  real::max,n
  !用矩阵(实型动态数组)将线性方程组表示出来!
  real,dimension(:,:),allocatable::a,b
  real,dimension(:),allocatable::x,m
  write(*,*)'请输入线性方程组的元数n:'
  read(*,*)t
  !给动态数组分配内存!
  allocate(a(t,t+1),m(t+1),x(t),b(t,t+1))
  write(*,*)'请依次按行输入线性方程组的增广矩阵a:'
  read(*,*)( b(i,:),i=1,t)
  a=b
  do k=1,t-1
    !选取每列的最大列元素!
    max=maxval(a(k:,k))!abs(a(k,k))
    imax=maxloc(a(k:,k),1) + k - 1    
    !将最大列元素所在的行与第K行进行交换!
    m(k:)    = a(k,k:)
    a(k,k:)    = a(imax,k:)
    a(imax,k:) = m(k:)
    !对方程组按X1,X2,...,Xt的顺序进行依次消元,将增广矩阵化为上三角矩阵!
    do i=k+1,t
      a(i,k+1:)=a(i,k+1:) - a(i,k)/a(k,k)*a(k,k+1:)
    end do
  end do
  !先算出Xt!
  x(t)=a(t,t+1)/a(t,t)
  !回带!
  do k=t-1,1,-1
    n= sum( a(k,k+1:t)*x(k+1:) )
    x(k)=(a(k,t+1)-n)/a(k,k)
  end do
  write(*,*)'线性方程组的解为:'
  write(*,*) x
  write(*,*) matmul( b(1:t,1:t) , reshape(x,[t,1]) )!//验证
  read(*,*)
end program S_equation

35

帖子

11

主题

0

精华

熟手

F 币
129 元
贡献
118 点
 楼主| 发表于 2017-2-23 18:04:31 | 显示全部楼层
fcode 发表于 2017-2-21 22:05
a(imax,j)=m(k,i)
改为
a(imax,j)=m(k,j)

谢谢大师兄的耐心,主要是你说的那个表控不熟悉,就不怎么敢用

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
发表于 2017-2-23 18:32:46 | 显示全部楼层
表控最简单了,尤其是 read,能表控坚决不用格式。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-4-19 08:08

Powered by Tencent X3.4

© 2013-2024 Tencent

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