Heat equation

This notebook showcases the solution of the heat equation for a few different cases

Heat equation in one dimention

The heatequation is given by: $\frac{\partial u}{\partial t} = \kappa \frac{\partial ^2 u}{\partial x^2}$

Case 1

With the given initial conditions: $u(0,t) = u(L,t) = 0$ and $u(x,o) = f(x)$, the solution to the heat equation is:

$u(x,t) = \sum_{n=1}^\infty b_n e^{(\frac{-\pi ^2 n^2}{L^2} \kappa t)} \cdot \sin{\frac{n\pi}{L}x}$

$b_n = \frac{2}{L} \int_0^L f(x) sin\frac{n\pi}{L}x$

Example 1

L = 50cm

f(x) = 100 deg celsius

$ b_n = \frac{200 (-1)^{n+1} + 200}{n\pi}$, which means that $b_n$ is 0 for even numbers and $400/n\pi$ for odd

$u(x,t) = \frac{400}{\pi} \sum_{n=0}^\infty \frac{1}{2n+1} e^{(\frac{-\pi ^2 (2n+1)^2}{(50)^2} \kappa t)} \cdot \sin{\frac{(2n+1)\pi}{50}x}$

With these initial conditions we can plot the function for any time t and position x

import numpy as np
import matplotlib.pyplot as plt

N = 100
XMIN = 0
XMAX = 50

YMIN = 0
YMAX = 1800

def fun(x, t, n):
    k=0.15
    pi2 = np.pi**2
    denom = 50**2
    kt = k*t
    
    return (400/np.pi) * ( (1/(2*n+1)) * np.exp(-(pi2 * (2*n+1)**2 * kt)/denom) * np.sin((2*n+1)*np.pi*x / 50) )


x = np.outer(np.linspace(XMIN, XMAX, N), np.ones(N))
y = np.outer(np.linspace(YMIN, YMAX, N), np.ones(N)).T
z = 0

for n in range(0,500):
    z += fun(x,y,n)

fig = plt.figure(figsize =(14, 9))
ax = plt.axes(projection ='3d')
ax.view_init(20,60)

ax.plot_surface(x, y, z)

<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x13837e0d0>

Example 2

L = 50cm

f(x) = 100 deg celsius

Now we can calculate the value for a specific position x and time t: For x = 25 and t = 1500 ($\kappa$ = 0.15 (copper))

$u(x,t) = \frac{400}{\pi} \sum_{n=0}^\infty \frac{1}{2n+1} e^{(\frac{-\pi ^2 (2n+1)^2}{(50)^2} \kappa (100))} \cdot \sin{\frac{(2n+1)\pi}{50}(1500)}$

import numpy as np
k=0.15
t=1500
x=25

pi2 = np.pi**2
denom = 50**2
kt = k*t

s=0
for n in range(0,1000):
    s += ( (1/(2*n+1)) * np.exp(-(pi2 * (2*n+1)**2 * kt)/denom) * np.sin((2*n+1)*np.pi*x / 50) )

u = 400/np.pi * s

print("Temperature at position " + str(x) + "cm and time " + str(t) + "s: " + str(u) + "(celcius)")

Temperature at position 25cm and time 1500s: 52.36282377966995(celcius)