The following simulation shows an example of the 1-D diffusion equation for the temperature in a 1-dimensional rod. The diffusion equation describes how the temperature T evolves with time t. For 1-D the diffusion equation is

This represents the temperature distribution in a plate of material. The temperature is shown using either a colour map or a graph of the intensity as a function of distance along the rod (or a combination of both). The temperature along the plate can be adjusted by clicking positions along the graph to adjust the intensity at that location. The temperature at the ends of the plate can be set using the sliders. The sides of the plate are fully insulated (no energy flows in or out). Click the start button to being the simulation.

Code Samples

Matlab

% Declare constants for the simulation
length = 0.2;           % Rod length (m)
plength = 100;          % Number of resolution elements
dx = length/plength;    % Position resolution (m)
D = 8.1e-5;             % Diffusion constant (SI, for aluminium)
dt = 1.0e-3;            % Time resultion (s)
Tl = 10.0;              % Temperature at the left boundary
Tr = 10.0;              % Temperature at the right boundary
outputRate = 1000;      % Number of steps between output cycles

% Set the initial temperature distribution
x = linspace(0.0, length, plength);
T = 80*sin(linspace(0.0, pi, plength)) + Tl; T = transpose(T);

% Create a matrix representation for our finite difference scheme
c = D*dt/dx/dx;
b = zeros(plength, 1); b(1) = Tl*c; b(end) = Tr*c;
A = diag(ones(plength, 1), 0)*(1.0 - 2.0*c) ...
    + diag(ones(plength-1, 1), 1)*c ...
    + diag(ones(plength-1, 1), -1)*c;

figure;
h = plot(x, T);
axis([0.0 length 0.0 100.0]);

% Evolve the system and show the animation
while ishandle(h)
    for ii=1:outputRate
        T = A*T + b;
    end
    set(h, 'ydata', T);
    pause(0.05);
    drawnow;
end

Python

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Declare constants for the simulation
length = 0.2;           # Rod length (m)
plength = 100;          # Number of resolution elements
dx = length/plength;    # Position resolution (m)
D = 8.1e-5;             # Diffusion constant (SI, for aluminium)
dt = 1.0e-3;            # Time resolution (s)
Tl = 10.0;              # Temperature at the left boundary
Tr = 10.0;              # Temperature at the right boundary
outputRate = 1000;      # Number of steps between output cycles

# Set the initial temperature distribution
x = np.linspace(0.0, length, plength);
T = 80*np.sin(np.linspace(0.0, np.pi, plength)) + Tl;
T = T.reshape(T.size,1)

# Create a matrix representation for our finite difference scheme
c = D*dt/dx/dx;
b = np.zeros((plength,1)); b[0] = Tl*c; b[-1] = Tr*c;
A = (np.diag(np.ones((plength,1))[:,0], 0)*(1.0 - 2.0*c)
    + np.diag(np.ones((plength-1,1))[:,0], 1)*c
    + np.diag(np.ones((plength-1,1))[:,0], -1)*c);

fig, ax = plt.subplots()
line, = ax.plot(x, T, lw=2)
ax.set_ylim(0.0, 100.0)
ax.set_xlim(0.0, length)

def animate(step_number):

    # Update the temperature distribution
    global T
    for i in xrange(0, outputRate):
        T = np.matmul(A, T) + b;

    line.set_data(x, T)
    return line,

ani = animation.FuncAnimation(fig, animate, xrange(1, 200),
    interval=25, blit=True, repeat=False)
#ani.save('output.mp4')
plt.show()

Cpp Gnuplot

#include <fstream>
#include <vector>
#include <cmath>

// Simulation constants
const double length = 0.2;            // Rod length (m)
const unsigned plength = 100;         // Number of resolution elements
const double dx = length/(plength-1); // Position resolution (m)
const double D = 8.1e-5;              // Diffusion constant (SI, for aluminium)
const double dt = 1.0e-3;             // Time resolution (s)
const double Tl = 10.0;               // Temperature at the left boundary
const double Tr = 10.0;               // Temperature at the right boundary
const unsigned outputRate = 1000;     // Number of steps between output cycles
const unsigned outputLength = 200;    // Length of simulation (# output files)
const double c = D*dt/dx/dx;          // Constant for the finite differences

int main(int argc, char** argv) {
  std::vector T(plength);

  // Set the initial temperature distribution
  for (unsigned k = 0; k < plength; ++k)
    T[k] = Tl + 80.0*sin((double)k/(plength-1)*M_PI);

  for (unsigned i = 0; i < outputLength; ++i) {
    for (unsigned j = 0; j < outputRate; ++j) {

      // Calculate the new temperature
      double tl = Tl;
      for (unsigned k = 0; k < plength-1; ++k) {
        double to = T[k];
        T[k] = c*tl + c*T[k+1] + (1.0 - 2.0*c)*T[k];
        tl = to;
      }
      T.back() = c*tl + c*Tr + (1.0 - 2.0*c)*T.back();
    }

    // Write the data out to a file
    char buf[128];
    snprintf(buf, 128, "output_%06d.dat", i);
    std::ofstream fp(buf);
    fp << "# x\tT\n";
    for (unsigned k = 0; k < plength; ++k) {
      fp << k*dx << '\t' << T[k] << '\n';
    }
  }

  return 0;
}

set terminal gif animate delay 10
set output 'output.gif'
set xrange [0.0:0.2]
set yrange [0.0:100.0]

set xlabel "Position (m)"
set ylabel "Temperature (T)"

do for [i=0:200-1] {
    plot sprintf("output_%06d.dat", i) with lines \
        title sprintf("Frame %d", i)
}

Temperature: 1-D Diffusion Equation

Code Samples