0
$\begingroup$

I am new to solving PDEs with finite difference methods. I implemented the code below to solve the heat equation following the explicit scheme, but when I plot the result I am suprised that the decay in time is very small. Could anyone tell me what I am doing wrong ?

Here is my code in Python. Below I did plot for four steps in time (2, 50, 100, 150) what the heat looks as a function in space :

from __future__ import division
import matplotlib.pyplot as plt
import numpy as np
import math

L =1.
N =200
J =200
k = 1./4.
dx = float(L)/float(N)
T = 0.5 * float(J)/k * dx * dx - 0.00001 #ensures that alpha <0.5
dt = float(T)/float(J)
alpha = k * dt / (dx*dx) # stability condition: alpha <0.5
print ("stability condition (alpha<0.5) --> alpha = ", alpha)
u = np.zeros((N,J))

def f(phi):
    return math.cos(np.pi*phi/float(L))

x = np.arange(0,L,dx)
t = np.arange(0,T,dt)

for n in range(0,N):
    u[n][0] = 6.*math.sin(np.pi*x[n]/float(L))

for j in range(0, J):
    u[0][j] = u[N-1][j] = 0.

for j in range(0, J-1):
    for n in range(1,N-1):
        u[n][j+1]  = (1.-2.*alpha)*u[n][j]+ alpha * ( u[n+1][j] +u[n-1][j])

u1 = []
u2 = []
u3 = []
u4 = []
x1 = []

for n in range(0,N):
    u1.append(u[n][2])
    u2.append(u[n][50])
    u3.append(u[n][100])
    u4.append(u[n][150])
    x1.append(x[n])

fig = plt.figure(0)
plt.plot(x1,u1)
plt.plot(x1,u2)
plt.plot(x1,u3)
plt.plot(x1,u4)
plt.show()
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.