Checking my Julia code solving the heat equation using the ultraspherical spectral method

I want to solve the usual heat equation:

\frac {\partial u}{\partial t}-\frac {\partial^{2}u}{\partial x^2}=0,

where u(x,t) represents the temperature at position x\in \left[0,1\right] and time t\in \left[0,1\right]. The thermal diffusivity is set to be 1. We assume homogeneous Dirichlet boundary conditions

u(0,t)=u(1,t)=0.

The initial condition is given by

u(x,0)=\text{exp}\left[-5\left(x-\frac {1}{2}\right)^{2}\right].

Step 1: Determine the splitting rank of the partial differential operator.

The partial differential operator (PDO) is given by

\mathcal {L}=\frac {\partial}{\partial t}-\frac {\partial^2}{\partial x^2},

whose splitting rank k=2 since we have the following decomposition:

\mathcal {L}=\left(\mathcal {D}\otimes \mathcal {J}\right)+\left(-\mathcal {J}\otimes \mathcal {D}^2\right),

where \mathcal {J} is the identity operator and \mathcal {D} is the first-order differential operator.

Step 2: Construct a generalized Sylvester matrix equation with an n_{t}\times n_{x} solution matrix.

Let’s ignore the boundary and initial conditions for now. Over the domain \left(x,t\right)\in\left[0,1\right]\times \left[0,1\right], we want to solve

\mathcal {L}u(x,t)=0.

By the above decomposition, the ultraspherical spectral method can be used to represent the above ordinary differential operators (ODOs).These operators can be truncated to derive the following generalized Sylvester matrix equation:

A_{1}XC^{T}_{1}+A_{2}XC^{T}_{2}=F,

where

A_{1}=\mathcal {P}_{n_{t}}\mathcal {D}_{1}\mathcal {P}^{T}_{n_{t}},\quad C_{1}=\mathcal {P}_{n_{x}}\mathcal {S}_{1}\mathcal {S}_{0}\mathcal {P}^{T}_{n_{x}},
A_{2}=\mathcal {P}_{n_{t}}\left(-\mathcal {S}_{0}\right)\mathcal {P}^{T}_{n_{t}},\quad C_{2}=\mathcal {P}_{n_{x}}\mathcal {D}_{2}\mathcal {P}^{T}_{n_{x}},

and F is the n_{t}\times n_{x} matrix of the bivariate Chebyshev expansion coefficients for the right-hand side, which is just the zero matrix. The operator \mathcal {S}_{0} converts a vector of Chebyshev coefficients to a vector of Ultraspherical(1) coefficients and, more generally, \mathcal {S}_{\lambda}, for \lambda\geq 1, converts a vector of Ultraspherical (\lambda) coefficients to a vector of Ultraspherical(\lambda+1) coefficients. For \lambda\geq 1, D_{\lambda} maps a vector of Chebyshev expansion coefficients to a vector of Ultraspherical(\lambda) expansion coefficients of the \lambda th derivative.

# Set the resolution levels for x and t
n_t=32;
n_x=32;

# Define the ODOs

#A1
D_1=Derivative();
D_1t=D_1:Chebyshev();
A_1=D_1t[1:n_t,1:n_t];

#C1
S_02=Conversion(Chebyshev(), Ultraspherical(2));
C_1=S_02[1:n_x,1:n_x];

#A2
S_01=Conversion(Chebyshev(), Ultraspherical(1));
A_2=-1*S_01[1:n_t,1:n_t];

#C2
D_2=D_1^2
D_2x=D_2:Chebyshev();
C_2=D_2x[1:n_x,1:n_x];

#F
F=zeros(n_t,n_x);

Step 3: Impose the boundary constraints.

Typically, the matrix equation A_{1}XC^{T}_{1}+A_{2}XC^{T}_{2}=F does not have a unique solution as the prescribed linear constraints \mathcal {B}_{x} and \mathcal{B}_{t} must also be incorporated. We have

u(x,t)\approx \sum_{i=0}^{n_{t}-1}\sum_{j=0}^{n_{x}-1}X_{ij}T_{i}\left(\psi(t)\right)T_{j}\left(\phi(x)\right),

where

\psi(t)=2t-1,\quad \phi(x)=2x-1.

On the left and right boundaries, we have the homogeneous Dirichlet boundary conditions u(0,t)=u(1,t)=0. By investigating the action of \mathcal {B}_{x} on \left\{T_{0}(\phi(x)),\cdots,T_{n_{x}-1}(\phi(x))\right\}, we can discretize any linear constraint of the form \mathcal {B}_{x}u(x,t)=\pmb{g}(t)=0 as

XB^{T}_{x}=G^{T},

where B_{x} is an 2\times n_{x} matrix and G is an 2\times n_{t} matrix where each row containing the first n_{t} Chebyshev coefficients of each component of \pmb{g}. In our case,

B_{x}=\begin{bmatrix}T_{0}(-1)&T_{1}(-1)&\cdots &T_{n_x -1}(-1)\\ T_{0}(1)&T_{1}(1)&\cdots &T_{n_x -1}(1) \end{bmatrix}=\begin{bmatrix}1&-1&1&\cdots \\ 1&1&1&\cdots \end{bmatrix},\quad G=\begin{bmatrix}0&0&\cdots &0\\ 0&0&\cdots &0 \end{bmatrix}.

Similarly, on the bottom boundary, we have the following initial condition:

u(x,0)=\text{exp}\left[-5\left(x-\frac {1}{2}\right)^{2}\right].

By investigating the action of \mathcal {B}_{t} on the basis \left\{T_{0}(\psi(t)),\cdots,T_{n_{t}-1}(\psi(t))\right\}, we can discretize \mathcal {B}_{t}u(x,t)=\pmb{k}(x) as

B_{t}X=K,

where B_{t} is an 1\times n_{t} matrix and K is an 1\times n_{x} matrix where each row containing the first n_{x} Chebyshev coefficients of each component of \pmb{k}. In our case,

B_{t}=\begin{bmatrix}T_{0}(-1)&T_{1}(-1)&\cdots &T_{n_t -1}(-1) \end{bmatrix}=\begin{bmatrix}1&-1&1&\cdots \end{bmatrix},\quad K\text{ contains the first $n_{x}$ Chebyshev coefficients of each component of $\pmb{k}$}.

To be compatible, we require

KB^{T}_{x}=B_{t}G^{T}.

By assumption, the prescribed linear constraints are linearly independent so the column ranks of B_{x} and B_{t} are 2 and 1, respectively. Without loss of generality, we further assume that the principal 2\times 2 and 1\times 1 submatrices of B_{x} and B_{t} are the identity matrices. Otherwise, permute the columns of B_{x} and B_{t}, and the corresponding rows/columns of X, so the principal 2\times 2 and 1\times 1 matrices \hat{B}_{x} and \hat{B}_{t} are invertible, then redefine as B_{x}\to \hat{B}^{-1}_{x}B_{x}, G\to \hat{B}^{-1}_{x}G, B_{t}\to \hat{B}^{-1}_{t}B_{t}, and K\to \hat{B}^{-1}_{t}K.

# Impose the boundary constraints

# Define a function that returns the initial condition at the bottom side
function η(x)
    exp(-5 * (x - 1/2)^2)
end

# Define B_t and K
domain = 0..1;
cheb_space = Chebyshev(domain);
cheb_η = Fun(t -> η(t), cheb_space);
B_t=zeros(1,n_t);
B_t[1,:]=[(-1)^(i-1) for i in 1:n_t];
K=zeros(1,n_x);
K[1,:]=cheb_η.coefficients[1:n_x];

# Define B_x and G
B_x=ones(2,n_x);
B_x[1,:]=[(-1)^(i-1) for i in 1:n_x];
G=zeros(2,n_t);
B_x=inv(B_x[1:2,1:2])*B_x;

Step 4: Solve the matrix equation with linear constraints.

Our approach is to use the linear constraints to remove degrees of freedom in X and thus obtain a generalized Sylvester matrix equation with a unique solution without constraints. We can modify the matrix equation to

\begin{align*} A_{1}XC^{T}_{1}+A_{2}XC^{T}_{2}-\left[\left(A_{1}\right)_{1:n_{t},1}B_{t}XC^{T}_{1}+\left(A_{2}\right)_{1:n_{t},1}B_{t}XC^{T}_{2}\right]=F-\left[\left(A_{1}\right)_{1:n_{t},1}KC^{T}_{1}+\left(A_{2}\right)_{1:n_{t},1}KC^{T}_{2}\right], \end{align*}

where we have used the constraint B_{t}X=K. Moreover, by rearranging, the left-hand side becomes

\begin{align*} \left[A_{1}-\left(A_{1}\right)_{1:n_{t},1}B_{t}\right]XC^{T}_{1}+\left[A_{2}-\left(A_{2}\right)_{1:n_{t},1}B_{t}\right]XC^{T}_{2}, \end{align*}

and since the principal matrix of B_{t} is the identity matrix, the first column of the matrix A_{j}-\left(A_{j}\right)_{1:n_{t},1}B_{t} for 1\leq j\leq 2 contains only zeros. Similarly, the condition XB^{T}_{x}=G^{T} can be used to further modify the matrix equation as follows:

\begin{align*} &\left[A_{1}-\left(A_{1}\right)_{1:n_{t},1}B_{t}\right]X\left[C_{1}-\left(C_{1}\right)_{1:n_{x},1:2}B_{x}\right]^{T}+\left[A_{2}-\left(A_{2}\right)_{1:n_{t},1}B_{t}\right]X\left[C_{2}-\left(C_{2}\right)_{1:n_{x},1:2}B_{x}\right]^{T}\\ &=F-\left[\left(A_{1}\right)_{1:n_{t},1}KC^{T}_{1}+\left(A_{2}\right)_{1:n_{t},1}KC^{T}_{2}\right]-\left[\left(A_{1}-\left(A_{1}\right)_{1:n_{t},1}B_{t}\right)G^{T}\left(C_{1}\right)^{T}_{1:n_{x},1:2}+\left(A_{2}-\left(A_{2}\right)_{1:n_{t},1}B_{t}\right)G^{T}\left(C_{2}\right)^{T}_{1:n_{x},1:2}\right], \end{align*}

so that the first two rows of the matrix \left[C_{j}-\left(C_{j}\right)_{1:n_{x},1:2}B_{x}\right]^{T} for 1\leq j\leq 2 contain only zeros.

Now, the first column of A_{j}-\left(A_{j}\right)_{1:n_{t},1}B_{t} and the first two rows of \left[C_{j}-\left(C_{j}\right)_{1:n_{x},1:2}B_{x}\right]^{T} are zeros, the matrix equation is independent of the first two columns and the first row of X. The right-hand side will be truncated accordingly by removing the last row and the last two columns. Therefore, the matrix equation can be reduced and then solved, obtaining a matrix X_{22}\in \mathbb {C}^{\left(n_{t}-1\right)\times \left(n_{x}-2\right)}, where

X=\begin{pmatrix} X_{11} & X_{12}\\ X_{21} & X_{22} \end{pmatrix},\quad X_{11}\in \mathbb {C}^{1\times 2},\quad X_{12}\in \mathbb {C}^{1\times (n_{x}-2)},\quad X_{21}\in \mathbb {C}^{(n_{t}-1)\times 2}.

Once we have computed X_{22}, we can recover X by using the linear constraints. For instance, since B_{t}X=K, X_{12}=K_{2}-B^{(2)}_{t}X_{22}, where K=\left[K_{1},K_{2}\right] with K_{1}\in \mathbb {C}^{1\times 2} and K_{2}\in \mathbb {C}^{1\times (n_{x}-2)}, and B_{t}=\left[B^{(1)}_{t},B^{(2)}_{t}\right] with B^{(2)}_{t}\in \mathbb {C}^{1\times (n_{t}-1)}. Furthermore, since XB^{T}_{x}=G^{T} and the 2\times 2 principal submatrix of B_{x} is the identity matrix, we have X_{21}=G^{T}_{2}-X_{22}\left(B^{(2)}_{x}\right)^{T}, where G=\left[G_{1},G_{2}\right] with G_{1}\in \mathbb {C}^{2\times 1} and G_{2}\in \mathbb {C}^{2\times (n_{t}-1)}, and B_{x}=\left[B^{(1)}_{x},B^{(2)}_{x}\right] with B^{(2)}_{x}\in \mathbb {C}^{2\times (n_{x}-2)}. Lastly, we can recover X_{11} using either of the two formulas

X_{11}=K_{1}-B^{(2)}_{t}X_{21},\quad X_{11}=G^{T}_{1}-X_{12}\left(B^{(2)}_{x}\right)^{T}.
# Update the ODOs and F
Ā_1=A_1-A_1[:,1]*B_t;
Ã_1=Ā_1[1:end-1,2:end];

C̄_1=C_1-C_1[:,1:2]*B_x;
C̃_1=C̄_1[1:end-2,3:end];

Ā_2=A_2-A_2[:,1]*B_t;
Ã_2=Ā_2[1:end-1,2:end];

C̄_2=C_2-C_2[:,1:2]*B_x;
C̃_2=C̄_2[1:end-2,3:end];

F̄=F-(A_1[:,1]*K*transpose(C_1)+A_2[:,1]*K*transpose(C_2))-((A_1-A_1[:,1]*B_t)*transpose(G)*transpose(C_1[:,1:2])
    +(A_2-A_2[:,1]*B_t)*transpose(G)*transpose(C_2[:,1:2]));
F̃=F̄[1:end-1,1:end-2];

Step 5: Solve a generalized Sylvester matrix equation.

We are left with a standard generalized Sylvester matrix equation of the form

\tilde{A}_{1}X_{22}\tilde{C}^{T}_{1}+\tilde{A}_{2}X_{22}\tilde{C}^{T}_{2}=\tilde{F}.

We can expand the matrix equation into an (n_{t}-1)(n_{x}-2)\times (n_{t}-1)(n_{x}-2) linear system

\left[\tilde{C}_{1}\otimes \tilde{A}_{1}+\tilde{C}_{2}\otimes \tilde{A}_{2}\right]\text{vec}\left(X_{22}\right)=\text{vec}\left(\tilde{F}\right),

where \otimes denotes the Kronecker product operator for matrices and vec$\left(C\right)$ denotes the vectorization of the matrix C formed by stacking the columns of C into a single column vector. Since we are using the ultraspherical spectral method, the matrices \tilde{A}_{j} and \tilde{C}_{j} are almost banded and hence, the matrix \sum_{j=1}^{2}\left(\tilde{C}_{j}\otimes \tilde{A}_{j}\right) is also almost banded.

# Solve the linear system

A=kron(sparse(C̃_1),sparse(Ã_1))+kron(sparse(C̃_2),sparse(Ã_2));
b=vec(F̃);
X̃_22=A\b;
X_22=reshape(X̃_22,n_t-1,n_x-2);
K_2=K[:,3:end];
K_1=K[:,1:2];
G_2=G[:,2:end];
G_1=G[:,1];
B̃_t=B_t[:,2:end];
B̃_x=B_x[:,3:end];
X_12=K_2-B̃_t*X_22;
X_21=transpose(G_2)-X_22*transpose(B̃_x);
X_11=K_1-B̃_t*X_21;
top_row = hcat(X_11, X_12);
bottom_rows = hcat(X_21, X_22);
X=vcat(top_row, bottom_rows);

Step 6: Form the solution.

Recall that

u(x,t)\approx \sum_{i=0}^{n_{t}-1}\sum_{j=0}^{n_{x}-1}X_{ij}T_{i}\left(\psi(t)\right)T_{j}\left(\phi(x)\right)=\sum_{i=0}^{n_{t}-1}\sum_{j=0}^{n_{x}-1}X_{ij}\cos\left(i\cos^{-1}\left(\psi(t)\right)\right)\cos\left(j\cos^{-1}\left(\phi(x)\right)\right),

where

\psi(t)=2t-1,\quad \phi(x)=2x-1.
# Form the solution

# Define the function ϕ(x)
function phi(x)
    2*x - 1
end

# Define the function ψ(t)
function psi(t)
    2*t - 1
end

x_range = range(0, 1, length = n_x);
t_range = range(0, 1, length = n_t);
T_i = zeros(n_t, n_t);
v=zeros(1,n_t);
for i in 1:n_t
    t_i = t_range[i]
    for j in 1:n_t
        v[j]=1
        T_i[i, j] = ChebyshevT(v)(psi(t_i))
        v=zeros(1,n_t)
    end
end
v=zeros(1,n_x);
T_j=zeros(n_x,n_x);
for j in 1:n_x
    x_j = x_range[j]
    for i in 1:n_x
        v[i]=1
        T_j[i,j]=ChebyshevT(v)(phi(x_j))
        v=zeros(1,n_x)
    end
end
u=T_i*(X*T_j);

The MATLAB code is straightforward, which is just

% Define the domain
dom = [0, 1, 0, 1];

% Create a chebfun2 of the initial condition
eta = chebfun(@(x) exp(-5*(x-1/2).^2), [0, 1]);

% Define the heat equation PDE
N = chebop2(@(u) diff(u, 1, 1) - diff(u, 2, 2), dom);

% Set homogeneous Dirichlet boundary conditions
N.lbc = 0;
N.rbc = 0;
N.dbc = eta;

% Solve the PDE
u = N\0;

With the Julia code, I get the following contour plot:


With the MATLAB code (which is the correct one), I get the following contour plot:

One thing I know for sure for my Julia code is that the initial condition I use does not match the boundary conditions, which I do not know how to deal with. However, this is not a problem for the MATLAB. But other than this issue, they’re still different contour plots.

I feel like I follow every step of the ultraspherical spectral method when I generate the Julia code, but I may made some mistakes that I did not realize. Any suggestions or comments are appreciated!

I haven’t looked at it in detail, but some of the conversion and derivative operators seem to be defined on the default domain, which is -1..1, instead of the one that you use in this problem, which is 0..1. Could you look into something like the example in Solving PDEs · ApproxFun.jl and change the domains as appropriate, and check if this helps?

Since the initial condition is defined on the original domain, which is [0,1], I use the this domain to get the Chebyshev coefficients. Then I define the function \phi(x)=2x-1 to change the domain from [0,1] to [-1,1].

Shouldn’t you be using the Chebyshev points? Because it looks like you are using equally spaced points

i.e.

x_range = (cos.(range(0, pi, length = n_x))-1)/2

I have

u(x,t)\approx \sum_{i=0}^{n_{t}-1}\sum_{j=0}^{n_{x}-1}X_{ij}T_{i}\left(\psi(t)\right)T_{j}\left(\phi(x)\right)=\sum_{i=0}^{n_{t}-1}\sum_{j=0}^{n_{x}-1}X_{ij}\cos\left(i\cos^{-1}\left(\psi(t)\right)\right)\cos\left(j\cos^{-1}\left(\phi(x)\right)\right),

where

\psi(t)=2t-1,\quad \phi(x)=2x-1.

This should work for any points for x,t between 0 and 1.

Yes, but all the advantages of using chebyshev polynomials come only when you use the correct set of nodes.

How are you computing the coefficients in the sum?

I have updated more details regarding this. Hope it’s helpful. For the bivariate Chebyshev series of u(x,t), the Chebyshev coefficients, which are the entries of X are computed by solving a generalized Sylvester matrix equation.

“Since the initial condition is defined on the original domain, which is [0,1], I use the this domain to get the Chebyshev coefficients. Then I define the function ϕ(x)=2x−1 to change the domain from [0,1] to [−1,1].”

But shouldn’t you then change the diffusion coefficient to 4 to get the same time behaviour ?

You could just follow this example:

Or there’s an example of Convection-Diffusion which could be modified here:

1 Like

Thanks! I saw both examples before, but I just want to see where my Julia code goes wrong. Theoretically, both MATLAB and Julia codes should give me the same result.

Why does the diffusion coefficient becomes 4? But I do feel like it has something to do with it. I tried 2, and I think it somehow works.

you do the change of variable x'=2x-1, so {\partial \over \partial x} = 2 {\partial \over \partial x'} and {\partial^2 \over \partial x^2} = 4 {\partial^2 \over \partial x'^2} . So you need to either change the diffusion coefficient or change the time variable when comparing the solution.

1 Like

I see. Since I also have the change of variable t'=2t-1, therefore the diffusion coefficient becomes two.

1 Like

yep, I did not see that you already changed time by a factor 2.