Heat equation

Heat equation#

Introduction#

Wait a few seconds until Python interaction ready! is shown in the top.

After that follow the instructions in the Grasple Exercise.

On this webpage the heat equation

\[ \frac{\partial u}{\partial t}=a^2\frac{\partial ^2u}{\partial x^2} \]

for \(x\in (0,L),\ t>0\) is considered, where \(a^2\) is is the thermal conductivity of a bar.

The initial and boundary conditions are

\[\begin{split} \begin{align*} u(x,0) &= f(x), \\ u(0,t) &= \phi(t), \\ u(L,t) &= \eta. \end{align*} \end{split}\]

This webpage constructs a (truncated) Fourier series of the exact solution \(u(x,t)\) of this problem using sympy, and then uses numpy and pyplot to display the obtained exact values.

The next cell imports a function that performs this.

import micropip
await micropip.install("ipympl")
import matplotlib
if not hasattr(matplotlib.rcParams, '_get'):
    matplotlib.rcParams._get = matplotlib.rcParams.get
%matplotlib widget
import sympy as sp
t,x = sp.symbols('t,x')
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from IPython.display import Markdown, clear_output, HTML
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
# import warnings
# warnings.filterwarnings('ignore',message='overflow encountered in exp')
# warnings.filterwarnings('ignore',message='overflow encountered in multiply')
# warnings.filterwarnings('ignore',message='invalid value encountered in multiply')
    
def Heat_equation(a,phi,eta,f,t_end,L,t,x):

    anew = a/L
    
    display(Markdown('# The animation is being rendered, please be patient.'))

    # symbols, functions and accuracy initialisation
    n = sp.symbols('n')
    g = sp.Function('g')(x)
    h = sp.Function('h')(x,t)
    v = sp.Function('v')(x)
    V = sp.Function('V')(x,t)
    N = 10

    # find function V that satisfies nonhomogeneous time-dependent BC's
    ode_v = sp.diff(v,x,x)
    sol = sp.dsolve(ode_v)
    constants = sp.solve([sol.rhs.subs(x,0)-phi,sol.rhs.subs(x,1)-eta],sp.symbols('C1 C2'))
    V = sol.rhs.subs(constants)

    # define new initial condition and right hand side
    g = f-V.subs(t,0)
    h = -V.diff(t)+anew**2*V.diff(x).diff(x)

    # find the Fourier sine expansion of h(x,t) with respect to x
    hn = sp.Function('hn')(t,n)
    hn = sp.integrate(h*sp.sin(n*sp.pi*x),(x,0,1))*2

    # find the Fourier sine expansion of g(x) with respect to x
    gn = sp.Function('gn')(n)
    gn = sp.integrate(g*sp.sin(n*sp.pi*x),(x,0,1))*2

    # find the Fourier sine expansion of w(t,x) with respect to x and construct w
    s = sp.symbols('s')
    Ln = sp.Function('Ln')(n)
    Ln = (anew*n*sp.pi)**2
    w = sp.Function('w')(x,t)
    w = 0
    for nn in np.arange(1,N+1):
        Lnn = Ln.subs(n,nn)
        hnn = hn.subs(n,nn)
        gnn = gn.subs(n,nn)
        wnn = sp.integrate(hnn.subs(t,s)*sp.exp(Lnn*s),(s,0,t))
        wnn = gnn+wnn
        wnn = wnn*sp.exp(-Lnn*t)
        w = w + wnn*sp.sin(nn*sp.pi*x)

    # make finally u(x,t)
    u = sp.Function('u')(x,t)
    u = w + V

    # evaluate u(x,t) at several x and t and plot
    u_lam = sp.lambdify((x,t),u,modules=['numpy'])
    t_vals = np.linspace(0,t_end,int(t_end)*100+1)
    x_vals = np.linspace(0,1,201)
    tv, xv = np.meshgrid(t_vals, x_vals)
    uv = u_lam(xv,tv)
    if isinstance(f,int):
        uv[:,0] = f
    else:
        fnum = sp.lambdify(x,f,"numpy")
        uv[:,0] = fnum(x_vals)
    uv = np.where(np.isinf(uv),np.zeros_like(uv),uv)
    uv = np.where(np.isnan(uv),np.zeros_like(uv),uv)
    
    u_min = np.min(np.min(uv))
    u_max = np.max(np.max(uv))
    
    fig = plt.figure(2);
    ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, L), ylim=(0,0.1*L));
    ax.grid();
    ax.set_xlabel('$x$');
    ax.set_aspect('equal','box')
    ax.set_yticks([])
    
    # test plot, afterwards, make []
    x_surf = np.vstack((x_vals,x_vals))*L
    y_surf = 0.1*L*np.ones_like(x_surf)
    y_surf[0,:] = 0
    u_surf = np.vstack((uv[:,0],uv[:,0]))
    clev = np.arange(u_min,u_max,0.001)
    surf = plt.imshow(u_surf, vmin = u_min, vmax = u_max, cmap=plt.cm.coolwarm, origin='lower', 
           extent=[0, L, 0, 0.1*L])
    plt.colorbar(label='$u(x,t)$',location='bottom')
    time_template = '$t$ = %1.3fs';
    time_text = ax.text(0.5, 1.05, '', transform=ax.transAxes,horizontalalignment='center');
    
    # u_surf = np.vstack((uv[:,-1],uv[:,-1]))
    # surf.set(data=u_surf)
    
    def init():
        u_surf = np.vstack((uv[:,0],uv[:,0]))
        surf = plt.imshow(u_surf, vmin = u_min, vmax = u_max, cmap=plt.cm.coolwarm, origin='lower', 
           extent=[0, L, 0, 0.1*L])
        time_text.set_text('');

        return surf, time_text
        
    def animate(i):
        u_surf = np.vstack((uv[:,i],uv[:,i]))
        surf.set(data=u_surf)
        time_text.set_text(time_template % (t_vals[i]));
        return surf, time_text
        
    skip = 0
    ani = animation.FuncAnimation(fig, animate, np.arange(0, len(t_vals),1+skip),
                                  interval=1, blit=True, init_func=None,repeat=False);
    plt.close(fig);
    
    ani_jshtml = ani.to_jshtml(fps=15)
    ani_html = HTML(ani_jshtml)
    clear_output()
    display(ani_html)
    display(Markdown('# You can slow down the animation (if you want) by using the - button.'))
    

    pass

INPUT CELL#

In the cell below you should put the values of the parameters you see in Grasple.

Please note: Each attempt in Grasple contains different values for the parameters.

After you run the cell the results of the calculations are shown with the set parameter values.

Blow the figure are controls that you can use to manipulate the animation.

''' Set the given values/functions in the following lines '''
a = sp.sqrt(0.3)
phi = 4*sp.sin(8*t)
eta = 0
f = 0
t_end = 2
L = 1

''' The next line performs the calculations and shows the results. '''
Heat_equation(a,phi,eta,f,t_end,L,t,x)
Contents