[Fortran] 纯文本查看 复制代码
program main
implicit none
integer,parameter::N=11
integer t
real*8::A(N),B(N),C(N),D(N),H(N)
integer,parameter::Dx=10,Dt=1
!空间和时间的步长
real,parameter:: k=0.1
real,parameter:: lemad=real(k*Dt)/real(Dx**2)
!含水层参数
integer i,j
!声明
do i=2,N-1
A(i)=-lemad
B(i)=1+2*lemad
C(i)=-lemad
end do
write(*,20)lemad
20 format(/,5X,'lemad=',F9.6)
do i =1,N
H(i)=10
end do
write (*,*) H(2:N-1)
H(1)=10
H(N)=12
write(*,15)
15 format(5X,'TIME',5X,'HEAD')!赋初值和边界值
t = 0
do j=1,101
t=t+Dt
do i=2,N-1
D(i)=H(i)
D(2)=D(2)+H(1)*lemad
D(N-1)=D(N-1)+H(N)*lemad
end do
call chase (N,A,B,C,D,H)
write(*,30) t,H(i)
30 format(1X,I3,5X,11F9.6)
end do
read(*,*)t
end program
subroutine chase (N,A,B,C,D,H)
!追赶法计算三对角线方程组
implicit none
integer:: N
real*8::A(N),B(N),C(N),D(N),H(N)
real*8::L(N),U(N)
integer::i,j
L(2)=B(2)
U(2)=D(2)/L(2)
do i = 3,N-1
L(i) = B(i)-A(i)*C(i-1)/L(i-1)
U(i) =(D(i)-A(i)*U(i-1))/L(i)
end do
!回代求H
H(N-1)=U(N-1)
do j= 2,N-3
i=N-1-j
H(i)=U(i)-C(i)*H(i+1)/L(i)
end do
return
end subroutine chase