Matlab: Timestep stability in a 1D heat diffusion model

Question

I have a 1D heat diffusion code in Matlab which I was using on a timescale of 10s of years and I am now trying to use the same code to work on a scale of millions of years. Obviously if I keep my timestep the same this will take ages to calculate but if I increase my timestep I encounter numerical stability issues.

My questions are:

How should I approach this problem? What affects the maximum stable timestep? And how do I calculate this?

Many thanks,

Alex

close all
clear all

dx = 4;    % discretization step in m
dt = 0.0000001; % timestep in Myrs
h=1000;        % height of box in m
nx=h/dx+1;
model_lenth=1; %length of model in Myrs
nt=ceil(model_lenth/dt)+1;     % number of tsteps to reach end of model
kappa = 1e-6; % thermal diffusivity
x=0:dx:0+h;     % finite difference mesh
T=38+0.05.*x;  % initial T=Tm everywhere ...
time=zeros(1,nt);
t=0;
Tnew = zeros(1,nx);

%Lower sill
sill_1_thickness=18;
Sill_1_top_position=590;
Sill_1_top=ceil(Sill_1_top_position/dx);
Sill_1_bottom=ceil((Sill_1_top_position+sill_1_thickness)/dx);

%Upper sill
sill_2_thickness=10;
Sill_2_top_position=260;
Sill_2_top=ceil(Sill_2_top_position/dx);
Sill_2_bottom=ceil((Sill_2_top_position+sill_2_thickness)/dx);

%Temperature of dolerite intrusions
Tm=1300;

T(Sill_1_top:Sill_1_bottom)=Tm; %Apply temperature to intrusion 1

% unit conversion to SI:
secinmyr=24*3600*365*1000000;   % dt in sec
dt=dt*secinmyr;

%Plot initial conditions
figure(1), clf
f1 = figure(1); %Make full screen
set(f1,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); 
plot (T,x,'LineWidth',2)
xlabel('T [^oC]')
ylabel('x[m]')
axis([0 1310 0 1000])
title(' Initial Conditions')
set(gca,'YDir','reverse');

%Main calculation
for it=1:nt

  %Apply temperature to upper intrusion
  if it==10;  
      T(Sill_2_top:Sill_2_bottom)=Tm; 
  end

  for i = 2:nx-1
      Tnew(i) = T(i) + kappa*dt*(T(i+1) - 2*T(i) + T(i-1))/dx/dx;
  end

  Tnew(1) = T(1);
  Tnew(nx) = T(nx);

  time(it) = t;

  T = Tnew; %Set old Temp to = new temp for next loop
  tmyears=(t/secinmyr);

  %Plot a figure which updates in the loop of temperature against depth
  figure(2), clf
  plot (T,x,'LineWidth',2)
  xlabel('T [^oC]')
  ylabel('x[m]')
  title([' Temperature against Depth after ',num2str(tmyears),' Myrs'])
  axis([0 1300 0 1000])
  set(gca,'YDir','reverse');%Reverse y axis

  %Make full screen     
  f2 = figure(2); 
  set(f2,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); 

  drawnow

  t=t+dt;
end

Answer 1

The stability condition for an explicit scheme like FTCS is governed by $r = K dt/dx^2 < 1/2$ or $dt < dx^2/(2K)$ where K is your coefficient of diffusion. This is required in order to make the sign of the 4th order derivative leading truncation error term be negative.

If you do not want to be limited by timestep I suggest using an implicit scheme (albeit at a higher of computational cost than an explicit scheme). This can be achieved simply by using backward Euler for the diffusion term instead of forward Euler. Another option is Crank-Nicholson which is also implicit.

Answer 2

@Isopycnal Oscillation is totally correct in that the maximum stable step is limited in an explicit scheme. Just for reference this is usually referred to as the discrete Fourier number or just Fourier number and can be looked up for different boundary conditions. also the following may help you for the derivation of the Implicit or Crank-Nicholson scheme and mentions stability Finite-Difference Approximations to the Heat Equation by Gerald W. Recktenwald.

Sorry I don't have the rep yet to add comments