CODEWITHSUNDEEP
pythonphysicspde
jcm
jcm

Reputation: 205

FTCS Algorithm for the heat equation

I am attempting to implement the FTCS algorithm for the 1 dimensional heat equation in Python. 

import numpy as np

L = 1 #Length of rod in x direction
k = 0.3 #Thermal conductivity of rod
tmax = 5 #how many seconds
nx = 100 #number of spacial steps
nt = 100; #number of time steps
xi = np.linspace(0,L,nx) 
ti = np.linspace(0,tmax,nt)
dx = L/(nx-1)
dt = tmax/(nt-1)
r = k*dt/(dx)**2
r2 = 1-2*r

u=np.zeros((nt,nx))
#IC
phi = 100;
for x in range(0,nx):    
    u[0][x] = phi

#BC
for t in range(0,nt):    
    u[t][0] = 0;
    u[t][nx-1] = 0
#FTCS Algorithm

for t in range(0,nt-1): #timestep
    for x in range(1,nx-2):
        u[t+1][x] = r*(u[t][x-1]+ u[t][x+1]) + r2*(u[t][x])

However I am not getting plausible values for u[t][x] = u(x,t) considering my initial conditions of 100. i.e. they blow up and give me stupid values like '4.11052068e+221' is there some bad programming practise I am partaking in which is destroying the algorithm? Or have I just implemented the algorithm incorrectly?

EDIT: I think I have worked out that it is because the algorithm is stable if and only if r < 1/2. The numbers just blow up because my r is about 2.5 or whatever, however if anyone can see any other errors let me know!!

Upvotes: 0

Views: 3163

Answers (1)

Mark Hannel
Mark Hannel

Reputation: 795

Original Question

The stability of the FTCS scheme hinges on the size of the constant r. If r<1/2, then rounding errors introduced at each step will exponentially decay. If r>1/2, then those rounding errors will exponentially increase. (As you've alluded to in your edit).

Small-ish Errors

  1. dx = L/nx and dt = tmax/nt. (If this is confusing to you, imagine the case of nx = 2 and L = 1... then dx = 0.5).
  2. L and tmax should be floats. (ie - L = 1.0).
  3. You are not updating the temperature u at the second to last position on the rod. Make sure the for loop over position extends to nx-1 and NOT nx-2.

Things to consider

  1. You might consider using real values/units for time, length, thermal conductivity and initial temperature.

  2. When you apply the boundary conditions (under the comment #BC), only the time t=0 has values of u which are non-zero, so you do not need to iterate through all the time steps.

  3. While it is necessary for you to iterate through time (hence the time for loop), you can actually vectorize the spatial derivative. You could replace your for loop with the following: for t in range(0, nt-1): u[t+1, 1:nx-1] = r*(u[t,0:nx-2] + u[t,2:]) + r2*(u[t,1:nx-1])

Upvotes: 1

Related Questions