子程序Wik运行失败，这个地方i一直在改变，但是x没有赋值
本帖最后由 huhelong 于 2023-12-21 10:46 编辑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为增广矩阵
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
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天的时间就足够学会了。 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))
谢谢已解决
页:
[1]