We are studying the 1D heat equation on an interval $[0, x_0] \subset \mathbb{R}$:
$$ \frac{\partial T}{\partial t}(x, t) = \chi \frac{\partial^2 T}{\partial x^2} (x,t) $$with initial conditions:
$$ T(x, t=0) = \left\{ \begin{array}{l l} 2 \, T_0 \frac{x}{x_0} & \text{if $0 \leq x \leq \frac{x_0}{2}$}\\ 2 \, T_0 \frac{x_0 - x}{x_0} & \text{if $\frac{x_0}{2} \leq x \leq x_0 $} \end{array} \right. $$There are two terms in this equation:
To discretize them, we are using a finite difference technique. To derive the discretization, we start with the well-known Taylor expansion of a function $f$ at point $a$ up to the second order in $h$, a variable that is taken as "small":
$$ f(a + h) = f(a) + h \, f'(a) + \frac{h^2}{2} f''(a) + o(h^2) $$Replacing $h$ by $-h$ yields the following:
$$ f(a - h) = f(a) - h \, f'(a) + \frac{h^2}{2} f''(a) + o(h^2) $$Summing the last two relations and reorganizing terms yields:
$$ f''(a) \approx \frac{1}{h^2} \left( f(a - h) - 2 f(a) + f(a + h) \right) $$Using terms up to the the first order yields the following discretization for a first order derivative in the neighborhood of a point $a$:
$$ f'(a) \approx \frac{1}{h} \left( f(a + h) - f(a) \right) $$Time and space are both discretized. This means that the temperature distribution we are solving for are indexed by integer indices for the space and the time variable.
We denote by $X^n$ the vector of unknowns at a given timestep (we exclude $T_1$ and $T_N$, as they are equal to zero):
$$ X^n = \begin{pmatrix} T_2^n \\ T_3^n \\ \vdots \\ T_{N-2}^n \\ T_{N-1}^n \end{pmatrix} $$Discretizing the original equation with this notation we can obtain multiple systems of equations to solve. We'll derive the most general case, using a parameter $\theta$ that acts as a weight for which values of time to consider in the spatial derivative. We obtain the following:
$$ \frac{X^{n+1} - X^n}{\delta t} = \frac{1 - \theta}{(\delta x)^2} M X^n + \frac{\theta}{(\delta x)^2} M X^{n+1} $$Where $\delta t$ and $\delta x$ are the mesh steps in time and space and $M$ a tridiagonal $N - 2 \times N -2$ size matrix with values $(1, -2, 1)$ on the tridiagonal.
Reordering the terms in the above equation yields the following:
$$ \left( I_n - \lambda \theta M \right) X^{n+1} = \left(I_n + \lambda (1 - \theta) M \right) X^n $$This last equation is a linear matrix system that can be solved with standard tools of linear algebra that ship with Numpy or Scipy. This is what we're going to do below.
We start by activating the pylab mode for this notebook.
%pylab inline
rcParams["figure.figsize"] = (10, 6)
Populating the interactive namespace from numpy and matplotlib
We then specify our space mesh.
meshsize = 101
x = linspace(0, 2 * pi, num=meshsize)
We then construct the tridiagonal matrix $M$ using the following snippet from StackOverflow.
import scipy.sparse
def get_matrix_M(meshsize):
diag_rows = np.array([[1] * meshsize,
[-2] * meshsize,
[1] * meshsize])
positions = [-1, 0, 1]
return scipy.sparse.spdiags(diag_rows, positions, meshsize, meshsize)
We verify that this works correctly below:
M = get_matrix_M(x.shape[0])
print M.shape
M.todense()
(101, 101)
matrix([[-2, 1, 0, ..., 0, 0, 0], [ 1, -2, 1, ..., 0, 0, 0], [ 0, 1, -2, ..., 0, 0, 0], ..., [ 0, 0, 0, ..., -2, 1, 0], [ 0, 0, 0, ..., 1, -2, 1], [ 0, 0, 0, ..., 0, 1, -2]])
We can now define the initial conditions as specified above.
def get_initial_conditions(x, T_0, x_0):
T = zeros(x.shape)
T[x<x_0/2] = 2 * T_0 * x[x<x_0/2]/x_0
T[x>=x_0/2] = 2 * T_0 * (1 - x[x>=x_0/2]/x_0)
return T
We plot this function to verify that we did code it well.
plot(x, get_initial_conditions(x, 1, x.max()))
[<matplotlib.lines.Line2D at 0x4a23e50>]
Now on to the next step: solving the linear system for which we have already constructed the matrix $M$. We proceed by building a class called Solution. The idea behind this is to initialize it with an initial state and calculate the matrix that solves the problem once, in the init procedure.
T_0 = get_initial_conditions(x, 1, x.max())
class Solution():
def __init__(self, x, T_0, lambda_value, theta):
self.x = x[1:-1]
self.T = matrix(T_0[1:-1]).T
self.t = 0
dx = x[1] - x[0]
self.dt = lambda_value * dx ** 2 # needed for updating the total time
right_hand_matrix = identity(self.T.shape[0]) + lambda_value * (1 - theta) * get_matrix_M(self.T.shape[0])
left_hand_matrix = identity(self.T.shape[0]) - lambda_value * theta * get_matrix_M(self.T.shape[0])
self.inv_matrix = dot(inv(left_hand_matrix), right_hand_matrix)
def solve_step(self):
self.T = self.inv_matrix * self.T
self.t += self.dt
def solve_n_steps(self, n):
for i in range(n):
self.solve_step()
sol = Solution(x, T_0, 0.5, 0.1)
sol.solve_n_steps(100)
plot(sol.x, array(sol.T), label='t = %.2f' % sol.t)
legend()
<matplotlib.legend.Legend at 0x4971d70>
Now that we have defined the problem and checked that our implementation behaves as expected, we're going to plot some solutions with different parameters to highlight different types of behaviors. To achieve this, we build a function that provides the solution at different timesteps and plots it.
def plot_iterations(theta, lambda_value, plot_legend=True):
sol = Solution(x, T_0, lambda_value, theta)
plot(sol.x, array(sol.T), label='t = %.2f' % sol.t)
sol.solve_n_steps(2)
plot(sol.x, array(sol.T), label='t = %.3f' % sol.t)
sol.solve_n_steps(8)
plot(sol.x, array(sol.T), label='t = %.2f' % sol.t)
sol.solve_n_steps(15)
plot(sol.x, array(sol.T), label='t = %.2f' % sol.t)
sol.solve_n_steps(25)
plot(sol.x, array(sol.T), label='t = %.2f' % sol.t)
sol.solve_n_steps(50)
plot(sol.x, array(sol.T), label='t = %.2f' % sol.t)
if plot_legend:
legend()
title(u"$\lambda=%.2f$ and $\\theta=%.1f$" %(lambda_value, theta))
We can now use that function to compare different parameter sets. Here, we compare the case $\lambda = 0.5$ to $\lambda = 1$ across the explicit Euler and the implicit Euler scheme ($\theta = 0$ and $\theta = 1$). As one can see, a value of $\lambda = 1$ is unstable with the explicit Euler scheme, while it is stable in the implicit Euler.
subplot(221)
plot_iterations(0, 0.5)
subplot(222)
plot_iterations(0, 1.0)
subplot(223)
plot_iterations(1, 0.5)
subplot(224)
plot_iterations(1, 1.0)
We can study the explicit Euler scheme more closely for parameters around $0.5$. We can really observe how the instability develops with this set of parameters.
subplot(221)
plot_iterations(0, 0.5, plot_legend=False)
subplot(222)
plot_iterations(0, 0.51, plot_legend=False)
subplot(223)
plot_iterations(0, 0.52, plot_legend=False)
subplot(224)
plot_iterations(0, 0.53, plot_legend=False)
Another thing to study is that stability does not necessarily mean convergence or good sampling. With the plots below, one sees that even though the scheme is stable (it doesn't explode), the solution can be not necessarily well discretized and thus exhibit "strange" behaviour, which disappears in the case of the implicit Euler scheme.
subplot(231)
plot_iterations(0.5, 20, plot_legend=False)
subplot(232)
plot_iterations(0.5, 60, plot_legend=False)
subplot(233)
plot_iterations(0.5, 100, plot_legend=False)
ylim(0, 1) #automatic scale somehow not good fit
subplot(234)
plot_iterations(1, 20, plot_legend=False)
subplot(235)
plot_iterations(1, 60, plot_legend=False)
subplot(236)
plot_iterations(1, 100, plot_legend=False)