# 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: 11 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: 219 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: 11 thanks a lot!

 January 27, 2015, 12:37 #5 New Member   Join Date: Jan 2015 Posts: 24 Rep Power: 11 Thanks a lot for the help! it was great!

 Tags crank nicolson, crank-nicolson, fortran code, ftn95, heat equation