I want to solve the usual heat equation:

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

The initial condition is given by

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

The partial differential operator (PDO) is given by

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

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

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:

where

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

where

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

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,

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

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

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,

To be compatible, we require

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

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

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:

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

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

```
# 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

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

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

where

```
# 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!