1D Diffusion Equation with FD+CN

diffusion1d.py

"""Solve the 1D diffusion equation using CN and finite differences."""
from time import sleep

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

# The total number of nodes
nnodes = 10
# The total number of times
ntimes = 500
# The time step
dt = 0.5
# The diffusion constant
D = 0.1
# The spatial mesh size
h = 1.0

G = nx.grid_graph(dim=[nnodes])
L = np.matrix(nx.laplacian(G))

# The rhs of the diffusion equation
rhs = -D*L/h**2

# Setting initial temperature
T = 60*np.matrix(np.ones((nnodes,ntimes)))
for i in range(nnodes/2):
    T[i,0] = 0;

# Setup the time propagator. In this case the rhs is time-independent so we
# can do this once.
ident = np.matrix(np.eye(nnodes,nnodes))
pmat = ident+(dt/2.0)*rhs
mmat = ident-(dt/2.0)*rhs
propagator = np.linalg.inv(mmat)*pmat

# Propagate
for i in range(ntimes-1):
    T[:,i+1] = propagator*T[:,i]

# To plot 1 time
# plot(T[:,10])

# To plot all times
# plot(T)