1DdiffusionTim

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()