1
$\begingroup$

I want to solve the heat equation $$u_t=au_{xx},\quad 0\leq t, \ 0\leq x\leq\pi$$ in Python using spectral methods. I set $u(0,x)=\sin(x)$ as initial value for the time and choose $a=2$. My Code:

import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint
from scipy.fftpack import diff as psdiff

def GenerateID(x):
    u = np.sin(x)
    return u

def hEq(u, t, L, a):
    uxx = psdiff(u, period=L,  order=2) 
    dudt = a*uxx
    return dudt

def hEq_solution(u0, t, L):
    sol = odeint(hEq, u0, t, args=(L,a,), mxstep=5000)
    return sol

if __name__ == "__main__":
    L = np.pi
    N = 1000
    x = np.linspace(0, L, N)
    a = 2
    u0 = GenerateID(x)
    Tfinal = 5.0
    t = np.linspace(0.0, Tfinal, 1000)
    sol = hEq_solution(u0, t, L)
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    sx, st = np.meshgrid(x, t)
    ax.plot_surface(sx, st, sol, cmap='jet')
    plt.xlabel('x')
    plt.ylabel('t')
    plt.show()

Which produces enter image description here This seems kind of right. Now my question: If I want to introduce the boundary condition $$u(t,0)=u(t,\pi)=0,$$ how would I do that? I know that considering all above conditions, the exact solution should be $$u(x,t)=\exp[-2t]\sin(x).$$

$\endgroup$
  • 3
    $\begingroup$ In a spectral method the boundary conditions are baked into the choice of the basis functions ($\sin(nx)$ in your case). Thus really the question is whether your fft is set up to do a sine transform or not which is now really a programming question. $\endgroup$ – Ian Mar 13 at 7:48
  • $\begingroup$ Please use the Math forum only for Math questions $\endgroup$ – David Mar 13 at 13:33
  • $\begingroup$ Well, it is a math question now. @Ian I don't understand how the basis would do that. Also, can I influence that by calculating the fft and change the wave numbers? $\endgroup$ – schoeni Mar 14 at 5:36
  • $\begingroup$ @schoeni The point is that you're hoping to decompose the initial condition into a sum of $\sin(nx)$ terms, since these are the eigenfunctions of the Laplacian with the boundary conditions you gave. This means a sine transform, not a Fourier transform. Then the heat equation is essentially an infinite diagonal system of ODEs once the initial data is represented in this basis. But whether or not your code does this depends entirely on what scipy.fftpack.diff does to the vector np.sin evaluated at np.linspace(0, np.pi,1000). $\endgroup$ – Ian Mar 24 at 2:08

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.