# Crank-Nicolson FTN95 code

 Register Blogs Members List Search Today's Posts Mark Forums Read January 5, 2015, 14:15 Crank-Nicolson FTN95 code #1 New Member   Join Date: Jan 2015 Posts: 24 Rep Power: 10 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 ************************************************** ***   January 6, 2015, 19:26 #2 Senior Member   Troy Snyder Join Date: Jul 2009 Location: Akron, OH Posts: 208 Rep Power: 18 I have a code that solves the 1d heat eqn. via Crank-Nicolson, I will compare tomorrow. nbarman likes this.   January 8, 2015, 10:52 #3 New Member   Join Date: Jan 2015 Posts: 24 Rep Power: 10 thanks a lot!   January 22, 2015, 10:14 #4 Member   Join Date: Dec 2012 Posts: 92 Rep Power: 12 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 Athos1387 likes this. Last edited by beer; January 23, 2015 at 07:23.   January 27, 2015, 12:37 #5 New Member   Join Date: Jan 2015 Posts: 24 Rep Power: 10 Thanks a lot for the help! it was great!  Tags crank nicolson, crank-nicolson, fortran code, ftn95, heat equation Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post holger_marschall OpenFOAM 242 March 7, 2013 12:30 cwang5 OpenFOAM Programming & Development 1 May 30, 2011 04:47 MechE OpenFOAM 28 May 16, 2011 11:02 Zdravko Stojanovic Main CFD Forum 2 July 19, 2010 10:11 morteza08 Main CFD Forum 5 June 1, 2010 08:21

All times are GMT -4. The time now is 13:22.