Skip to content

Instantly share code, notes, and snippets.

@thider
Created May 18, 2011 18:13
Show Gist options
  • Save thider/979154 to your computer and use it in GitHub Desktop.
Save thider/979154 to your computer and use it in GitHub Desktop.
1DdiffusionTim
"""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 = np.matrix(np.eye(nnodes,nnodes))
D = .1*D
# this loop trys to set a few of the central D values to be different
#so we can simulate something being stuck in the center of the device
Drod = 2
for i in range(4):
D[i + nnodes/2, i + nnodes/2] = 1
D[i + nnodes/2, i + nnodes/2] = Drod*D[i + nnodes/2, i + nnodes/2]
print D
# 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[nnodes/2,i] = T[nnodes/2, i] + 1
T[:,i+1] = propagator*T[:,i]
# To plot 1 time
# plt.plot(T[:,100])
# plt.show()
# To plot all times
plt.plot(T)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment