CRANK-NICOLSON SCHEME TO SOLVE HEAT DFFUSION EQUATIONI

[watto8]

QUESTION: Heat diffusion equation is u_t= (D(u)u_x)_x.
Modify this program to investigate the following developments:
Allow for the diffusivity D(u) to change discontinuously, with initial data as u(x,0)= (1+x)(1-x)^2. boundary values u(+-1,t)=0.
Treat in detail the case D(u)=1 when x<1/2 and D(u)=1/2 otherwise.
=> This is my normal code:

 program crank_nicolson
implicit none
real, allocatable :: x(:),u(:),a(:),b(:),c(:),d(:)
real:: m,dx,dt,tmax
integer:: 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
ji= 20
m = dt/(2*dx**2)

allocate (x(0:ji),u(0:ji+1),a(ji),b(ji),c(ji),d(ji))
open(10,file='crank_nicolsan.m')

!initial condition
do j=1,ji
  x(j)= -1+j*dx
  u(j)=(1+x(j))*(1-x(j))**2 !x=-1+j*dx
end do
x(0)= -1
x(JI)= 1

!boundary condition
u(0)=0
u(ji+1)=0 

  do j = 1, ji ! a,b,c are the coefficients of c-n scheme and d is right part
    a(j) = -m
    b(j) = 1+2.0*m
    c(j) = -m
    d(j) = m*u(j+1)+(1-2.0*m)*u(j)+m*u(j-1) 
  end do
  call thomas(a,b,c,d)
  do j=1,ji-1
    u(j)=d(j)
  end do
 
!Print out the Approximate solution in matlab file 
write(10,*)  'ApproximateSolution =[',x(0),u(0)
do j =1, JI-1  
write(10,*) x(j),u(j)
 end do
write(10,*)x(JI),u(JI),']'
 write(10,*) "plot(ApproximateSolution(:,1),ApproximateSolution(:,2),'r')"
write(10,*) "xlabel('x'),ylabel('temperature'),legend('Approximate C-N Scheme')"
close(10) 
contains

subroutine thomas (a,b,c,d)

  real, intent(inout) :: a(:),b(:),c(:),d(:)
  integer j

  do j = 2,ji !combined decomposition and forward substitution
!	 a - sub-diagonal (means it is the diagonal below the main diagonal)
!	 b - the main diagonal
!	 c - sup-diagonal (means it is the diagonal above the main diagonal)
!	 d - right part
    
    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(ji) = d(ji)/b(ji)
  do j = ji-1,1,-1
    d(j) = (d(j)-c(j)*d(j+1))/b(j)
  end do

end subroutine thomas 

end program crank_nicolson

(THERE ARE NO ERRORS IN THE ABOVE CODE. JUST NEED TO MODIFY THE ABOVE CODE Allowing the diffusivity D(u) to change discontinuously, with initial data as u(x,0)= (1+x)(1-x)^2. boundary values u(+-1,t)=0.
Treat in detail the case D(u)=1 when x<1/2 and D(u)=1/2 otherwise

 

[FJacq]

Sorry but they are errors in your code. You have modified wrongly what I gave you. And they are also mistakes in my version...

Now it is time to think by yourself. First of all, find the relationship between the diffusivity and your variable m.

François Jacq

 

[watto8]

Francios Jacq, Its a different question but just need to need the previous code to solve this new questions. As fae as Im concerned, this code is working perfectly.

 

[watto8]

Treating the case
D(u)=1 when x<1/2
D(u)=1/2 When x>1/2 (otherwise)

I need to use if statement:

if (-1+j*dx >1/2) then
.....


end if ..
but not sure where to put and there might other codes I need to add to get discontinuity when x>1/2 at D(u)=1/2

and from get the values of u(j) values and plot the graphs. And I think the graphs will contains some solid curve and also some dashes line to show discontinuity.