子程序Wik运行失败，这个地方i一直在改变，但是x没有赋值

huhelong

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

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

Subroutine uptri(Wik,ft,F,N)!!A(N,N)系数矩阵,b(N)右向量，N方程维数，F方程的根
Implicit real*8(a-z)
Integer::c,d,N
real*8::Wik(N,N),ft(N),F(N)
F(N)=ft(N)/Wik(N,N)
  do c=N-1,1,-1
      F(c)=ft(c)
      do d=c+1,N
          F(c)=F(c)-Wik(c,d)*F(d)
      end do
      F(c)=F(c)/Wik(c,c)
  end do 
end Subroutine uptri
end Module m_gauss

Subroutine Wik_(G,M,W,x,y,z,H)
Implicit none
Integer::a=0,b=0,L=1,i,k,j,N=11
real*8::G(11),M(11),W(11,11),x(11),y(11),z(11),H,pi=acos(-1.0)
Do i=1,N
   x(i)=(L/2)*(1-cos(((i-1)/(N-1))*pi))
END Do
Do k=1,N
   z(k)=(L/2)*(1-cos(((k-1)/(N-1))*pi))
END Do
Do j=1,N
   y(j)=(L/2)*(1-cos(((j-1)/(N-1))*pi))
END Do 
 DO
100 a=a+1
Do  
b=1+b
if (a==b)then
   i=a
   k=b
    G(i)=product(x(i)-y,y/=x(i)) 
    W(i,k)=1/G(i)
    write(*,*)"W(",i,"",k,")=",W(i,k)
end if
if(a/=b.and.b<=11)then
    i=a
    k=b
    G(i)=product(x(i)-y,y/=x(i))  
    M(k)=product(z(k)-y,y/=z(k))
    H=1/(x(i)-x(k))
    W(i,k)=H*G(i)/M(k)
write(*,*)"W(",i,"",k,")=",W(i,k)
end if
if (b>11)then
    b=0
    Goto 100
    end if
end do
if(a==12)then
exit
end if
end do
end Subroutine Wik_   

program main_mik
use m_gauss
Implicit none
Integer,parameter::N=11
integer::c,d
real*8::Wik(11,11),W(11,11),ft(11),F(11),G(11),M(11),x(11),y(11),z(11),H
call Wik_(G,M,W,x,y,z,H)
Wik=W
ft=(/1,1,1,1,1,1,1,1,1,1,1/)
call lineq(Wik,ft,F,N)
open (unit=12,file='fout.txt')
write (12,101)F
101format(T5,'计算结果',/,T4,'F=',4(/F12.8))
end program main_mik

necrohan

Integer::a=0,b=0,L=1,i,k,j,N=11
Do i=1,N
   x(i)=(L/2)*(1-cos(((i-1)/(N-1))*pi))
END Do
这里 L 是整数1，L/2=0，所以计算结果一直是0。如果要不是0，L需要声明为实数。
程序里对y，z的循环也一样，如果y，z和x相同，可以直接赋值 y=x; z=x;

另外，大哥，求你先学一点fortran的基础知识，先找本书看看，我看了你的其他提问，都特别基础，有1天的时间就足够学会了。

huhelong

necrohan 发表于 2023-12-21 13:04
Integer::a=0,b=0,L=1,i,k,j,N=11
Do i=1,N
   x(i)=(L/2)*(1-cos(((i-1)/(N-1))*pi))

谢谢已解决