Problem 9.12

Conclusion

When running the implicit FTCS method and comparing to fig 6.7 and 6.8, found that when tau smaller than tsigma, the method is stable as with explicit FTCS fig 6.7. But when tau > tsigma, the explicit FTCS method (fig 6.8) showed it is unstable, however the implicit method showed that it is still stable even though tau > tsigma as can be seen by the plot below. So implicit FTCS is always stable.

When running C-N method, it was also stable when tau> tsigma. The plots show C-N is similar to implicit FTCS method.

part a

OUTPUT

» nma_problem_9_12

function to solve problem 9.12

Nasser Abbasi

Enter N, number of grid points (suggest 61):61

Enter tau, the time step (suggest 1e-4 or 1.5e-4)1e-4

Enter nstep, the number of steps to run (suggest 300):300

Enter method 1=implicit FTCS, 2=Crank-Niclson:1

»

» nma_problem_9_12

function to solve problem 9.12

Nasser Abbasi

Enter N, number of grid points (suggest 61):61

Enter tau, the time step (suggest 1e-4 or 1.5e-4)1.5e-4

Enter nstep, the number of steps to run (suggest 300):300

Enter method 1=implicit FTCS, 2=Crank-Niclson:1

»

part b

» nma_problem_9_12

function to solve problem 9.12

Nasser Abbasi

Enter N, number of grid points (suggest 61):61

Enter tau, the time step (suggest 1e-4 or 1.5e-4)1e-4

Enter nstep, the number of steps to run (suggest 300):300

Enter method 1=implicit FTCS, 2=Crank-Niclson:2

»

» nma_problem_9_12

function to solve problem 9.12

Nasser Abbasi

Enter N, number of grid points (suggest 61):61

Enter tau, the time step (suggest 1e-4 or 1.5e-4)1.5e-4

Enter nstep, the number of steps to run (suggest 300):300

Enter method 1=implicit FTCS, 2=Crank-Niclson:2

»

code:

function nma_problem_9_12()

%

%function to solve problem 9.12

%Nasser Abbasi

%

% This is the same as dftc except we use the implicit

% Crank-Nicolson or the implicit FTCS method for

% solving the diffusion problem

%

clear all; help nma_problem_9_12;

L=1;

kappa=1; % diffusion coeff.

N=input('Enter N, number of grid points (suggest 61):');

h=L/(N-1);

tau=input('Enter tau, the time step (suggest 1e-4 or 1.5e-4)');

tsigma=h^2/(2*kappa);

coeff1= tau/(2*tsigma); % for FTCS

coeff2= tau/(4*tsigma); % for C-N

nstep=input('Enter nstep, the number of steps to run (suggest 300):');

method = input('Enter method 1=implicit FTCS, 2=Crank-Niclson:');

% make the D matrix based on Boundary conditions of Derchlit

D=makeD(N);

tt=zeros(N,1);

tt(round(N/2))=1/h ; %initial cond is delta at center

xplot = (0:N-1)*h - L/2; % xscale

iplot = 1; %counter for plots

nplots =50; % number of snapshots

plot_step=nstep/nplots;

% calculate the B matrix and use it in the loop

I=eye(N);

if(method ==1) % implict FTCS

B= inv(I-(coeff1*D));

else

B= inv(I - (coeff2*D)) * (I + (coeff2*D));

end

for istep=1:nstep

tt= B*tt; % this is the implicit method

if(rem(istep,plot_step)<1)

ttplot(:,iplot)=tt(:);

tplot(iplot)=istep*tau;

iplot = iplot+1;

end

end

figure;

mesh(tplot,xplot,ttplot);

xlabel('Time');

ylabel('x');

zlabel('T(x,t)');

if(method ==1)

title('Diffusion of delta spike using the implict FTCS method');

else

title('Diffusion of delta spike using the C-N method');

end

pause(1);

figure;

contourLevels=0:0.5:10;

contourLabels=0:5;

cs=contour(tplot,xplot,ttplot,contourLevels);

clabel(cs,contourLabels);

xlabel('Time');

ylabel('x');

title('Temperature contour plot');

%%%%%%%%%%%%%%%%%%%%%%%%%

%

%

%%%%%%%%%%%%%%%%%%%%%%%%%

function D=makeD(N)

D=zeros(N);

for(i=2:N-1)

D(i,(i-1))= 1;

D(i,(i)) = -2;

D(i,(i+1))= 1;

end