Crank-Nicolson FTN95 code
 

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

http://www.axgig.com/images/57869547764844014047.jpg

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
************************************************** ***

I have a code that solves the 1d heat eqn. via Crank-Nicolson, I will compare tomorrow.

thanks a lot!

Do you still need help? I found the code quite interesting using the Thomas algorithm, so I've taken a look.

I think your coefficients are not quite right. Think about the equation again, separate the time steps n+1 and n and then build the coefficient matrix. At the Neuman boundary condition I used the shortcut that zero gradient means zero flux, so you don't need a dummy point outside the domain to reconstruct the gradient. Other than that your code seems all right and after fixing the coefficients the solution looks ok I guess.

Code:

program HW36u
implicit none
integer, parameter:: mag = 21
real , dimension(mag) :: u , uold , x , d , b , dprime
real , dimension(2:mag) :: a
real , dimension(mag-1) :: cprime , c
real , dimension( mag , mag ) :: 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 = 1. / (mag-1)
r = dt / (dx**2)
u0 = 0 !boundary condition

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

!! Build coefficient matrix for steady state solution
!constructing the diagonals of the coeficents matrix A
do i=1,mag
b(i) = 2*r
end do
b(mag) = r
do i=1,mag-1
c(i) = -r
end do
c(1) = 0
do i=2,mag
a(i) = -r
end do
!constructing the coeficents of knwon u, matrix K
do i=1,mag
do j=1,mag
if ( i == j ) then
K(i,j) = -b(i)
else if ( j-i == 1) then
K(i,j) = -c(i)
else if ( i-j == 1 ) then
K(i,j) = -a(i)
else
K(i,j) = 0
end if
end do
end do

!! Now add transient term on both sides
do i=1,mag
b(i) = b(i) + 2
K(i,i) = K(i,i) + 2
end do

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

!Thomas Algorithm
do m = 1,n

uold = u

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

!constructing c'
do i=1,mag-1
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,mag
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=mag,1,-1
if ( i == mag ) then
u(i) = dprime(i)
else
u(i) = dprime(i) - ( cprime(i)*u(i+1) )
end if
end do
end do


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

!output data to a file
open (1, file='hw3-6.txt')
do i=1,mag
write(1,100) x(i), ' ', u(i)
end do
100 format (F6.3,A3,F6.3)
close(1)

end program HW36u



Solution:
x u
0.000 0.000
0.050 0.021
0.100 0.042
0.150 0.063
0.200 0.086
0.250 0.110
0.300 0.135
0.350 0.162
0.400 0.190
0.450 0.219
0.500 0.248
0.550 0.279
0.600 0.308
0.650 0.337
0.700 0.365
0.750 0.390
0.800 0.412
0.850 0.431
0.900 0.446
0.950 0.456
1.000 0.461

Thanks a lot for the help! it was great!