求助：一维非定常热传导方程BTCS(隐式格式)的求解Fortran程序
我是个程序新手，但是现在急需要求解一维非定常热传导方程的隐式格式，需要永fortran或者matlab的程序，盼望好心的高手帮帮忙啊！不生感激！ program btcs
c...Performs backward-time central-space(IMPLICIT EULER)
parameter (nmax=100,delta=.000001)
double precision lambda,u(-2:nmax+2),h(1:nmax),u0(1:nmax+1)
double precision a(1:nmax), b(1:nmax), c(1:nmax)
double precision d(1:nmax), e, f
open(unit=9,file='btcs.out')
c...Read initial data samples.Samples evenly spaced.
c...Data assumed periodic.
open(unit=8,file='nb.dat',status='old')
read(8,*) n, lambda, tfinal
if(n.gt.nmax) then
write(9,*) '****Too many data points****'
close(unit=8)
close(unit=9)
stop
endif
if(n.lt.2) then
write(9,*) '****Too few data points****'
close(unit=8)
close(unit=9)
stop
endif
if(lambda.lt.0.01) then
write(9,*) '****Lambda small or negative****'
close(unit=8)
close(unit=9)
stop
endif
i=1
read(8,*,err=1000,end=1000) xmin, u(1)
do 10, i=2,n
read(8,*,err=1000,end=1000) dummy, u(i)
10 continue
i=n+1
read(8,*,err=1000,end=1000) xmax, u(n+1)
if(abs(u(n+1)-u(1)).gt..0001) then
write(9,*) '****Data not periodic****'
close(unit=8)
close(unit=9)
stop
endif
if(xmax.le.xmin+.0001) then
write(9,*) '****Bad x-axis****'
close(unit=8)
close(unit=9)
stop
endif
u(0) = u(n)
u(-1) = u(n-1)
u(-2) = u(n-2)
u(n+2) = u(2)
do 15, i=1,n+1
u0(i) = u(i)
15 continue
delta_x=(xmax-xmin)/real(n)
delta_t=lambda*delta_x
itert=nint(tfinal/delta_t)
write(9,*) 'Final time requested: ', tfinal
tfinal = real(itert)*delta_t
write(9,*) 'Actual final time: ', tfinal
write(9,*) 'delta_t = ', delta_t
write(9,*) 'delta_x = ', delta_x
write(9,*) 'lambda = ', lambda
do 500, it=1,itert
c...Form ``periodic tridiagonal'' system of equations
c...BTCS
do 20, i= 1,n
c...Lower sub-diagonal
a(i) = -.5*lambda
c...Diagonal
b(i) = 1.0
c...Upper sub-diagonal
c(i) = .5 * lambda
c...Right hand side
d(i) = u(i)
c...Upper right element
e = -.5* lambda
c...Lower left element
f = .5 * lambda
20 continue
c...BTBS
c do 20, i= 1,n
c...Lower sub-diagonal
c a(i) = -lambda
c...Diagonal
c b(i) = 1.0 + lambda
c...Upper sub-diagonal
c c(i) = 0.
c...Right hand side
c d(i) = u(i)
c...Upper right element
c e = - lambda
c...Lower left element
c f = 0.
c20 continue
call ptrid(n,a,b,c,d,e,f,h)
do 150, i=1,n
u(i) = h(i)
150 continue
u(0) = u(n)
u(-1) = u(n-1)
u(-2) = u(n-2)
u(n+1) = u(1)
u(n+2) = u(2)
500continue
write(9,*)
write(9,1050)
do 800, i=1,n+1
write(9,1100) i, u0(i), u(i), abs(u(i)-u0(i))
800continue
close(unit=8)
close(unit=9)
C...write simple file for plotting
open(unit=10,file='btcs.plt')
do 900, i=1,n+1
write(10,1150) -1 + real(2*i-2)/real(n), u(i),u0(i)
900continue
close(unit=10)
stop
1000 write(9,*) '****Error reading data point number ', i,'****'
close(unit=8)
close(unit=9)
stop
1050 format(4x,'N',11x,'INITIAL',12x,'BTCS',11x, 'DIFFERENCE')
1100 format(i5,5x,f13.8,5x,f13.8,5x,f13.8)
1150 format(f15.9,5x,e25.9,5x,f15.9)
end
subroutine ptrid(n,a,b,c,d,e,f,x)
integer n
double precision a(n),b(n),c(n),d(n),e,f,x(n)
double precision r(n),s(n)
r(1) = e
s(1) = f
c...Zero out the subdiagonal
do 10, i=2,n-2
c a(i) = a(i) - a(i)*b(i-1)/b(i-1)
b(i) = b(i) - a(i)*c(i-1)/b(i-1)
d(i) = d(i) - a(i)*d(i-1)/b(i-1)
r(i) = - a(i)*r(i-1)/b(i-1)
10 continue
c a(n-1) = a(n-1) - a(n-1)*b(n-2)/b(n-2)
b(n-1) = b(n-1) - a(n-1)*c(n-2)/b(n-2)
c(n-1) = c(n-1) - a(n-1)*r(n-2)/b(n-2)
d(n-1) = d(n-1) - a(n-1)*d(n-2)/b(n-2)
c a(n) = a(n) - a(n)*b(n-1)/b(n-1)
b(n) = b(n) - a(n)*c(n-1)/b(n-1)
d(n) = d(n) - a(n)*d(n-1)/b(n-1)
if(abs(f).gt..0001) then
c...Chase the non-zero element from left to right across the
c...bottom row.
do 15, i=1,n-2
c s(i) =s(i) - s(i)*b(i)/b(i)
s(i+1) = - s(i)*c(i)/b(i)
b(n) =b(n) - s(i)*r(i)/b(i)
d(n) =d(n) - s(i)*d(i)/b(i)
15 continue
c s(n-1) =s(n-1) - s(n-1)*b(n-1)/b(n-1)
b(n) =b(n) - s(n-1)*c(n-1)/b(n-1)
d(n) =d(n) - s(n-1)*d(n-1)/b(n-1)
endif
c...Solve the upper-triangular system using
c...back-substitution.
x(n) = d(n)/b(n)
x(n-1) = (d(n-1)-c(n-1)*x(n))/b(n-1)
do 20, i=n-2,1,-1
x(i) = (d(i)-c(i)*x(i+1)-r(i)*x(n))/b(i)
20 continue
return
end 感谢您的帮助！:handshake
页:
[1]