[Athos1387]

Hi,
I must solve the question below using crank-nicolson method and Thomas algorithm by writing a code in fortran.



I've written a code for FTN95 as below. It works without a problem and gives me the answers, the problem is that the answers are wrong. I've solved it with FTCS method and analytically,and I know what the right answers are. I've checked the code 100 times and I don't understand what's wrong with it. can anyone pleeeeeas help me?!!!

************************************************** **
program HW36u
implicit none
real , dimension(19) :: u , uold , x , d , b , dprime
real , dimension(2:19) :: a
real , dimension(18) :: cprime , c
real , dimension( 19 , 19 ) :: K
real :: t , dt , dx , r , u0 , u20
integer :: j , i , n , m
write (*,*)'type your desirable time'
read (*,*) t
write (*,*)'type your desirable time step'
read (*,*) dt
n = t / dt
dx = 0.05
r = dt / (dx**2)
u0 = 0 !boundary condition

!Initial Condition
do i=1,19
x(i) = 0 + (i)*dx
if ( x(i) < 0.75 ) then
uold(i)= 0
else if ( x(i) > 0.75 ) then
uold(i) = 1
else
uold(i)= 0.5
end if
end do

!constructing the diagonals of the coeficents matrix A
do i=1,19
b(i) = (-2) * (r+1)
end do
do i=1,18
c(i) = r
end do
do i=2,19
a(i) = r
end do


!constructing the coeficents of knwon u, matrix K
do i=1,19
do j=1,19
if ( i == j ) then
K(i,j) = (2) * (r-1)
else if ( j-i == 1 .or. i-j == 1 ) then
K(i,j) = (-1) * r
else
K(i,j) = 0
end if
end do
end do

!computing the value of heat eq. for the desired time

!Thomas Algorithm
do m = 1,n

!constructing the knwon matrix d
d = matmul (K , uold)

!constructing c'
do i=1,18
if (i==1) then
cprime(i) = c(i) / b(i)
else
cprime(i) = c(i) / ( b(i) - cprime(i-1)*a(i) )
end if
end do

!constructing d'
do i=1,19
if (i==1) then
dprime(i) = d(i) / b(i)
else
dprime(i) = ( d(i) - dprime(i-1)*a(i) ) / ( b(i) - cprime(i-1)*a(i) )
end if
end do

!calculating u
do i=19,1,-1
if ( i == 19 ) then
u(i) = dprime(i)
else
u(i) = dprime(i) - ( cprime(i)*u(i+1) )
end if
end do

uold = u

end do

!Neuman boundary condition
u20 = ( (real(4)/3)*u(19) ) - ( (real(1)/3)*u(18) )

!output on screen
write (*,*) 'The data will be outputed to a file'

!output data to a file
open (1, file='hw3-6.txt' , status='new')
write(1,*) u0
do i=1,19
write(1,*) u(i)
end do
write(1,*) u20
close(1)

end program HW36u
************************************************** ***