求助Gauss消元法求解线性方程组的数值不稳定问题

Yjc

各位老师好，我在学习数值算法的时候写了一个用Gauss消元法求解线性方程组的Fortran代码，并重复了徐士良老师《Fortran常用算法程序集》中的例子（见下图）；这段代码编译过程正常，但求解出的X1和X2与书中给出的答案相差很大。其中一个原因我认为是选主元方式的差异，我的程序中采用“列选主元”，而书中采用“全选主元”。很多教材中提到列主元也能在很大程度上保证数值稳定性，所以我想原因可能出在代码的其他部分，希望各位老师帮忙看看这段代码在算法上有没有可以改进的地方，感谢各位老师不吝指教！
（经过测试，发现是回代这一部分循环的代码出现问题，因为我用直接列式算出的结果是正确的，但我实在没有发现回代这个循环的问题）

2.f90 (1.33 KB, 下载次数: 4)

2024-3-28 12:29 上传

点击文件名下载附件
源代码

Yjc

Transpose 发表于 2024-3-28 13:55
[mw_shl_code=fortran,true]   !ans=0 这个ans的位置错误
    X(4) = B(4)/A(4,4)
    X(3) = (B(3)-A(3,4 ...

感谢老师指正，我搞错了ans赋值语句的位置，现在程序运行正常了。
还想麻烦老师帮忙看看这个程序在数值算法和语法规范等方面有没有可以优化的地方，非常感谢！

Transpose

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

   !ans=0 这个ans的位置错误 
    X(4) = B(4)/A(4,4)
    X(3) = (B(3)-A(3,4)*X(4))/A(3,3)
    X(2) = (B(2)-A(2,4)*X(4)-A(2,3)*X(3))/a(2,2)

    do j=3,1,-1
        ans = 0 !应该在这里
        do i=j+1,4,1
            ans = ans + A(j,i)*X(i)
        end do
        X(j) = (B(j)-ans)/A(j,j)
    end do