- #1

jkthejetplane

- 29

- 4

- Homework Statement
- I was tasked with putting a few methods for solving PDEs into matlab. I started with the given code that solved using FTCS method. Below is was i have modified it to. I understand i need to have it run FTCS for the first step. I am not sure if the way i have written the equations into the code makes sense but i have gotten the Richardson method to produce results. There is an error with the Dufort method though

"Array indices must be positive integers or

logical values.

Error in problem3 (line 45)

tt(2:(N-1)) = tt(2:(N-1)) +

2*coeff(tt(3:N) + tt(1:(N-2)) -

tt(2:(N-1)))/(1+ 2*coeff);

"

- Relevant Equations
- FTCS equation is the same form is copied below from my text

FTCS scheme

My Matlab code:

%Problem 3

%Solve diffusion problem using Richardson scheme or DuFort-Frankel scheme

clear all;

scheme = menu('Choose method of solving diffusion equation:', 'Richardson', 'DuFort-Frankel');

tau = input('Enter time step: ');

N = input('Enter the number of grid points: ');

L = 1.; % The system extends from x=-L/2 to x=L/2

h = L/(N-1); % Grid size

k = 1.; % Diffusion coefficient

coeff = k*tau/h^2;

if( coeff < 0.5 )

disp('Solution is expected to be stable');

else

disp('WARNING: Solution is expected to be unstable');

end

%* Set initial and boundary conditions.

tt = zeros(N,1); % Initialize temperature to zero at all points

tt(round(N/2)) = 1/h; % Initial cond. is delta function in center

% The boundary conditions are tt(1) = tt(N) = 0

%* Set up loop and plot variables.

xplot = (0:N-1)*h - L/2; % Record the x scale for plots

iplot = 1; % Counter used to count plots

nstep = 300; % Maximum number of iterations

nplots = 50; % Number of snapshots (plots) to take

plot_step = nstep/nplots; % Number of time steps between plots

%Run First steps using FTCS

tt(2:(N-1)) = tt(2:(N-1)) + coeff*(tt(3:N) + tt(1:(N-2)) - 2*tt(2:(N-1)))

%Run selected scheme

if (scheme==1) %richardson

for istep=2:nstep

tt(2:(N-1)) = tt(2:(N-1)) + 2*coeff*(tt(3:N) + tt(1:(N-2)) - 2*tt(2:(N-1)));

if( rem(istep,plot_step) < 1 ) % Every plot_step steps

ttplot(:,iplot) = tt(:); % record tt(i) for plotting

tplot(iplot) = istep*tau; % Record time for plots

iplot = iplot+1;

end

end

else (scheme==2)%dufort

for istep=2:nstep

tt(2:(N-1)) = tt(2:(N-1)) + 2*coeff(tt(3:N) + tt(1:(N-2)) - tt(2:(N-1)))/(1+ 2*coeff);

if( rem(istep,plot_step) < 1 ) % Every plot_step steps

ttplot(:,iplot) = tt(:); % record tt(i) for plotting

tplot(iplot) = istep*tau; % Record time for plots

iplot = iplot+1;

end

end

end

%* Plot temperature versus x and t as wire-mesh and contour plots.

figure(1); clf;

mesh(tplot,xplot,ttplot); % Wire-mesh surface plot

xlabel('Time'); ylabel('x'); zlabel('T(x,t)');

title('Diffusion of a delta spike');

pause(1);

figure(2); clf;

contourLevels = 0:0.5:10; contourLabels = 0:5;

cs = contour(tplot,xplot,ttplot,contourLevels); % Contour plot

clabel(cs,contourLabels); % Add labels to selected contour levels

xlabel('Time'); ylabel('x');

title('Temperature contour plot');