Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

CRANK-NICOLSON SCHEME TO SOLVE HEAT DFFUSION EQUATIONI

Status
Not open for further replies.

watto8

Programmer
Feb 3, 2014
29
GB
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:

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
 
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
 
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.
 
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.
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top