Transpose 发表于 2023-11-30 16:45
41行

谢谢，已解决

fcode 发表于 2023-12-1 08:25
1. 建议每个 module 下面都写 implicit none
2. 建议不要写任何   Implicit real*8(a-z) ，也就是说 implic ...

谢谢您的建议

1. 建议每个 module 下面都写 implicit none
2. 建议不要写任何   Implicit real*8(a-z) ，也就是说 implicit 后面只写 none
3. 对contains在module下的过程（函数和子程序），不写 implicit（它会继承module的implicit none）。
4. 没有contains在module和program下面的过程（外部过程），写 implicit none

41行

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

Module m_gauss
Contains
Subroutine lineq(A,b,x,N)!!解线性方程组 高斯列主元消去法Ax=b
  Implicit real*8(a-z)
  Integer::i,k,N
  Integer::id_max !主元素符号
  real*8::A(N,N),b(N),x(N)!!A(N,N)系数矩阵,b(N)右向量，N方程维数，x方程的根
  real*8::Aup(N,N),bup(N)
  real*8::Ab(N,N+1)!Ab为增广矩阵[Ab]
  real*8::vtemp1(N+1),vtemp2(N+1)
  Ab(1:N,1:N)=A
  Ab(:,N+1)=b
  do k=1,N-1
      elmax=dabs(Ab(k,k))
      id_max=k
!!!查找主元素（不是赋值最大元素给elmax，而是为了找出最大元素对应的标号）
      do i=k+1,n
       if (dabs(Ab(i,k))>elmax) then
           elmax=Ab(i,k)
           id_max=i
       end if
      end do
!!!完成查找最大元素，查找完成以后与第k行交换，其他不变
   vtemp1=Ab(k,:)
   vtemp2=Ab(id_max,:)
   Ab(k,:)=vtemp2
   Ab(id_max,:)=vtemp1
!!!交换两行元素，交换完成后按照消元法进行
    do i=k+1,N
       temp=Ab(i,k)/Ab(k,k)
       Ab(i,:)=Ab(i,:)-temp*Ab(k,:)
    end do 
  end do
  Aup(:,:)=Ab(1:N,1:N)
  bup(:)=Ab(:,N+1)
 call uptri(Aup,bup,x,n)!!!调用上三角方程组的回带方法
end Subroutine lineq

Subroutine uptri(A,b,x,N)!!A(N,N)系数矩阵,b(N)右向量，N方程维数，x方程的根
Implicit real*8(a-z)
Integer::i,j,k
real*8::A(N,N),b(N),x(N)
x(N)=b(N)/A(N,N)
  do i=n-1,1,-1
      x(i)=b(i)
      do j=i+1,N
          x(i)=x(i)-a(i,j)*x(j)
      end do
      x(i)=x(i)/A(i,i)
  end do 
end Subroutine uptri
end Module m_gauss

Module m_newton
    Contains
    Subroutine solve()!!计算非线性方程组，需要用到线性方程组的解法，这里用主元消元法；需要准备函数文件和偏导数文件
    use m_gauss
    Implicit real*8(a-z)
    Integer::I,itmax=50 !itmax最大允许迭代次数
    Integer::N=2  !!N方程维数  
    real*8::x(2),f(2),dx(2)
    real*8::df(2,2)!!df偏导数矩阵
    open(unit=11,file='result.txt')
    write(11,101)
101 format(/,T6,'牛顿法计算非线性方程组迭代序列',/)    
    x=(/2d0,2d0/)
    tol=1d-8
    do i=1,itmax
        call func(f,x)
        call jac(df,x)
        call lineq(df,-f,dx,N)
        x=x+dx
        write(11,102)i,x
102  format(I5,2F16.10)       
     dx2=dsqrt(dx(1)**2+dx(2)**2)!判断计算精度，当满足误差容限时退出循环。deqrt是 平方根
     if(dx2<tol)exit
    end do
    end Subroutine solve

    Subroutine func(f,x)!f方程函数，x自变量.
    Implicit real*8(a-z)
    real*8::x(2),f(2)
    f(1)=6*x(1)**3+x(1)*x(2)-3*x(2)**3-4
    f(2)=x(1)**2-18*x(1)*x(2)**2+16*x(2)**3+1
    end Subroutine func

    Subroutine jac(df,x)
    Implicit real*8(a-z)
    real*8::x(2),df(2,2)!!偏导数矩阵
    df(1,1)=18*x(1)**2+x(2)
    df(2,1)=2*x(1)-18*x(2)**2
    df(1,2)=x(1)-9*x(2)**2
    df(2,2)=-36*x(1)*x(2)+48*x(2)**2
    end Subroutine jac
    end Module m_newton


Program main
use m_newton 
call solve()
end Program main