# Help implementing finite difference scheme for heat equation

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.

• def u(x, t): return u, while it happens to be valid python, is not at all meaningful here. Perhaps it might be easier to start with a beginner programming textbook of some sort? – Kirill Jul 5 '16 at 2:05
• Javier can you please post your final code – user22301 Nov 10 '16 at 8:56

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.

• Hey, just wanted to say thanks so much for your help. You were very helpful in getting my program to work :) – Javier Jul 12 '16 at 0:31
• @Javier no problem man! – spektr Jul 12 '16 at 0:37

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.

• Hi, thank you so much. Your comment has been very helpful. I am still having some issues though. First, may I define the solution as follows: $u = np.zeros((N,M))$? I think it defines an $NxM$ matrix with zeros for all entries, but I am not sure. Then my loop should just fill in the values as it iterates.... Also you said I can set the initial condition with $u[:,0] = x**2$. So can I similarly define the B.C. with $u[0,:] = 0$ and $u[N,:] = 0$? When I run it this way I get the error that for $u[N,:] = 0$ the "index 10 is out of bounds for axis 0 with size 10". How can I correct this? – Javier Jul 6 '16 at 19:46
• Hey, no worries man! Yeah, using numpy.zeros is one way to do it. What you wrote actually creates an array of size $N \times M$ so the indices go from 0 to $N-1$ and $M-1$ respectively. Therefore, indices $N$ or $M$ are out of bounds :) – Chris Jul 11 '16 at 12:42
• I actually got it to work. Thanks so much for your help :) – Javier Jul 12 '16 at 0:30