Homework & Coursework Questions

! The one-dimensional PDE for heat diffusion equation
! u_t=(D(u)u_x)_x + s where u(x,t) is the temperature,
! D(u) is the diffusivity and s(x,t) is a source term.
! Taking D(u)= 1 and s(x,t)=0 gives
! u_t= u_xx
! uniform one dimensional region |x|<1 for t>0
! uniform mesh size delta x=0.1
! Initial condition : u(x)= (1+x)(1-x)^2
 
Program crank_nicolson
Implicit none
Real, allocatable :: x(:),u(:,:),a(:),b(:),c(:),d(:)
Real:: m,dx,dt,Tmax
Integer:: n,j,NI,JI
Print*, 'Enter the total number of time steps'
read*, NI
Print*, 'Enter the final time'
read*, Tmax
dt=Tmax/NI !The size of timestep
Print*, 'This gives stepsize of time dt=',dt 
dx= 0.1 !delta x =0.1
allocate (x(0:JI),u(0:NI,0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI))
m = dt/(2*dx**2)
JI= 20
open(10,file='crank_nicolsan.m')
!Initial Condition
do j=1,JI-1
x(j)= -1+j*dx
u(0,j)=(1+x(j))*(1-x(j))**2 !x=-1+j*dx
end do
x(0)= -1
x(JI)= 1
!Boundary condition
do n=0,NI
u(n,0)=0 !b.c. at j=0
u(n,JI)=0 !b.c. at j=JI
end do
do n =0 , NI !loop
do j = 1, JI ! a,b,c,d are the coefficients of C-N scheme 
a(j) = -m
b(j) = 1+2*m
c(j) = -m
d(j) = m*u(n,j+1)+(1-2*m)*u(n,j)+m*u(n,j-1) 
end do
end do
 
do j=1, JI-1
do n=1, NI-1
u(n,j)=d(j)
end do
end do
call thomas(a,b,c,d,JI)
!Print out the Approximate solution in matlab file
write(10,*) 'ApproximateSolution =[',x(0),u(0,0)
do j =1, JI
write(10,*) x(j),u(n,j)
end do
write(10,*)x(JI),u(NI,JI),']'
 
end Program crank_nicolson
 
subroutine thomas (a,b,c,d,JI)
implicit none
real, intent(inout) :: a(*),b(*),c(*),d(*)
integer, intent(in) :: JI 
integer j
 
do j = 2,JI !combined decomposition and forward substitution
a(j) = a(j)/b(j-1)
b(j) = b(j)-a(j)*c(j-1)
d(j) = d(j)-a(j)*d(j-1)
end do 
 
!back substitution
d(j) = d(j)/b(j)
do j = JI-1,1,-1
d(j) = (d(j)-c(j)*d(j+1))/b(j)
end do
return
end subroutine thomas

allocate (x(0:NI),u(0:JI),a(0:JI),b(0:JI),c(0:JI),d(0:JI))