У меня есть код Fortran, который вычисляет вектор решения с помощью подпрограммы алгоритма Томаса.
Я хочу, чтобы вектор решения выполнялся в al oop в течение определенного количества времени.
Как мне вызвать эту подпрограмму в l oop?
Моя подпрограмма - это подпрограмма алгоритма Томаса.
Она возвращает вектор решения u, но я хочу, чтобы векторы использовались NN раз в аль oop. Таким образом, старый u становится новым u для использования в подпрограмме.
Как мне это сделать? Ниже я попробовал
program thomasalg2
implicit double precision(A-H,O-Z)
real*8, dimension(9,1) :: a,b,c,r,u,uold!the dimension is subject to change depending on the size of the new matrix
!real*8, dimension(9,50) :: W
real*8 :: pi
real*8 :: h,k,lm,l,T
integer :: i,j,al,NN,n
l = 1!right endpoint on the X-axis
n = 9 !number of rols/cols of the coefficient matrix with boundaries included
T = 0.5 !maximum number of the time variable
NN = 50!number of time steps
np = n
h = l/n
k = T/NN
al = 1.0D0 !alpha
pi = dacos(-1.0D0)
lm = (al**2)*(k/(h**2)) !lambda
do i = 1,n
r(i,1) = sin(pi*i*h) !this is W_0
end do
a(1,1) = 0.0D0
do i = 2,n
a(i,1) = -lm
end do
do i = 1,n
b(i,1) = 1 + (2*lm)
end do
c(9,1) = 0.0D0
do i = 1,n-1
c(i,1) = -lm
end do
!the 3 diagonals are stored in the 1st, 2nd, 3rd & 4th files respectively.
open(10, file = 'thom1.txt')
open(11, file = 'thom2.txt')
open(12, file = 'thom3.txt')
open(13, file = 'thom4.txt')
write(10,*)
do i = 1,n
write(10,*) a(i,1)
end do
write(11,*)
do i = 1,n
write(11,*) b(i,1)
end do
write(12,*)
do i = 1,n
write(12,*) c(i,1)
end do
write(13,*)
do i = 1,n
write(13,*) r(i,1)
end do
open(14, file = 'tridag2.txt')
write(14,*)
n = 9
do i = 1,n
write(14,*) a(i,1),b(i,1),c(i,1),r(i,1) !write the given vectors in the file in the form of a column vector
end do
call tridag(a,b,c,r,u,n)
!solve the given system and return the solution vector u
do i = 1,NN
call tridag(a,b,c,r,u,n)
!write(15,*) u
r = u
end do
open(15, file = 'tridag2u.txt')
write(15,*)
!write the solution vector in the form of a column vector
do i = 1,n
write(15,*) u(i,1)
end do
!print *, "Your data has been written in 'tridag2.txt'"
конец программы thomasalg2
подпрограмма tridag (a, b, c, r, u, n) неявной двойной точности (AH, OZ)
integer n, NMAX
real*8 a(n), b(n), c(n), r(n), u(n)
parameter (NMAX = 500)
integer j
real*8 bet, gam(NMAX)
if(b(1).eq.0.) stop "tridag: rewrite equations"
bet = b(1)
u(1)=r(1)/bet
do j = 2,n
gam(j) = c(j-1)/bet
bet = b(j)-a(j)*gam(j)
if (bet.eq.0.) stop "tridag failed"
u(j) = (r(j)-a(j)*u(j-1))/bet
end do
do j = n-1,1,-1
u(j) = u(j)-gam(j+1)*u(j+1)
end do
!print *, "The solution is", u
return
конец подпрограммы