Help implementing finite difference scheme for heat equation

Question

I am trying to solve the following problem via a finite difference approximation:

$u_t = k \, u_{xx}$, on $0 < x < L$ and $t > 0$;

$u(0,t) = u(L,t) = 0$;

$u(x,0) = f(x)$.

I take $u(x,0) = f(x) = x^2$ for my problem.

Programming is not my forte at all, so I am having some issues with the implementation. Here is my code:

 ## This program is to implement a Finite Difference method approximation
## to solve the Heat Equation, u_t = k * u_xx,
## in 1D w/out sources & on a finite interval 0 < x < L. The PDE
## is subject to B.C: u(0,t) = u(L,t) = 0,
## and the I.C: u(x,0) = f(x).
import numpy as np
import matplotlib.pyplot as plt

# definition of solution to u_t = k * u_xx
def u(x,t):
    return u

# definition of initial condition function  
def f(x):
    return x^2

# parameters    
L = 1
T = 10
N = 10
M = 10
s = 0.25

# uniform mesh
x_init = 0
x_end = L
dx = float(x_end - x_init) / N

# time discretization
t_init = 0
t_end = T
dt = float(t_end - t_init) / M

t = np.zeros(M+1)
t = np.arange(t_init, t_end, dt)

# Boundary Conditions
for m in xrange(0, M):
    t[m] = m * dt
    u(0, t[m]) = 0
    u(N, t[m]) = 0

# Initial Conditions
for j in xrange(0, N):
    x[j] = j * dx
    u(x[j], 0) = f(x[j])

# finite difference scheme
for j in xrange(1, N-1):
    u(x[j],t[m+1]) = u(x[j],t[m]) + s * ( u(x[j+1],t[m]) - 
    2 * u(x[j],t[m]) + u(x[j-1],t[m]) )

So in particular, I am wondering (1) if my function definitions are okay? (2) is my uniform mesh and time discretization okay? (3) am I selecting reasonable values for the parameters? (I know that the approximation is more accurate for smaller $\Delta x$ and $\Delta t$, with the disadvantage of increased computing time..)

(4) how are my loops (for the initial/boundary conditions)? I get the feeling I need a nested loop for the finite difference scheme.. One for loop for the values of j and another for loop going through each value m in time. Is this correct? ... Also, I keep getting the error: u(0, t[m]) = 0 "can't assign to function call." What is wrong there?

Thanks in advance for the help! As I mentioned, my programming skills are not very strong so I feel like a fish out of water doing this.

Answer 1

So your first big issue is one you bring up in point (4), where you say you get the error 'u(0, t[m]) = 0 "can't assign to function call"'. You're trying to store data by assigning data to a function. That won't work. The way you're implementing this code, you actually should make u(x,t) be represented by a 2-D array that can be accessed via u[m][j], where m is associated with the $m^{th}$ time value, and j is associated with the $j^{th}$ spatial value.

As you also note, you should have a nested for loop. Do an outside loop for time, and then an inner loop for looping through the spatial discretization.

Now assuming your value s makes sense for your problem, the other potential issue I see is your step size, dt. I am pretty sure this time stepsize won't be stable based on the value of dx since it seems you're just using an Explicit Euler time stepping scheme. Try picking some dt that is proportional to dx*dx, like dt=0.1*dx*dx, and see how it fairs for you.

Answer 2

Don't worry about programming skills, everyone's a beginner at some point. It's a good start actually. I'm just going to give some hints assuming that this is an exercise. Actually, choward's answer is complete, this is just a rephrasing of his suggestions in, I hope, a more beginner friendly language. So, if you do accept an answer it should be his ;)

1. Your definition of the solution does not make sense in Python. I would suggest storing it in a NxM numpy array. That way you can easily define the initial and boundary conditions with the slice notation (i.e. U[:,0] = x**2 where x is an 1-dim array, see numpy.linspace).

2. It seems that you are using the central difference scheme for the second order derivative in space. The definition is usually $$ u_{xx}(x_i) = \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} $$ where $h$ is the step size in space (your dx).

3. You are using a explicit time stepping scheme called explicit Euler or forward Euler which is conditionally stable. This means that your discretisation has to fulfil a certain condition in order for your solution to remain bounded or, in other words, that your time step has to be small enough. Another around this would be to use an implicit method such as backward Euler which is unconditionally stable. Note that, while this is fine for the heat equation which is intrinsically stable, it's maybe not the best choice for neutral or unstable problems since your numerical solution may suffer from over-stabilisation (a.k.a. numerical damping).

4. You have only implemented the inner loop. You would need an outer loop to advance in time as well.