Sylvester equation

• discretized Poisson equation:

\[ \mathbf{D}_{xx}\mathbf{U} + \mathbf{U}\mathbf{D}_{yy}' = \mathbf{F} \]

• Kronecker product, \(\otimes\):

A = [1 2; -2 0]
B = [ 1 10 100; -5 5 3 ]
A_kron_B = [ A[1,1]*B  A[1,2]*B;
             A[2,1]*B  A[2,2]*B ]
kron(A,B)

• Poisson as a linear system:

\[ [(\mathbf{I}_y \otimes \mathbf{D}_{xx}) + (\mathbf{D}_{yy} \otimes \mathbf{I}_x)]\text{vec}(\mathbf{U}) = \text{vec}(\mathbf{F})\\ \mathbf{Au} = \mathbf{b}. \]

• add boundary conditions to their precise values as in collocation strategy
• template: create the linear system, modify it for the boundary conditions, solve it using backslash, and reshape to get a grid function
  • condition number improved by rescaling by the largest of linear coefficients
• walk through poissonfd

ex:

f = (x,y) -> x^2 - y + 2;
m,n = 6,5
x,Dx,Dxx = FNC.diffmat2(m,[0,3]);
y,Dy,Dyy = FNC.diffmat2(n,[-1,1]);
unvec = u -> reshape(u,m+1,n+1);

F = [ f(x,y) for x in x, y in y ]
A = kron(I(n+1),Dxx) + kron(Dyy,I(m+1));
b = vec(F);

@show N = length(F);

spy(sparse(A),color=:blues,m=3,
    title="System matrix before boundary conditions")
isboundary = trues(m+1,n+1);
isboundary[2:m,2:n] .= false;
idx = vec(isboundary);

spy(sparse(isboundary),m=3,color=:darkblue,legend=:none,
    title="Boundary points",
    xaxis=("column index",[0,n+2]),yaxis=("row index",[0,m+2]) )

I_N = I(N);
A[idx,:] .= I_N[idx,:];     # Dirichlet conditions

spy(sparse(A),color=:blues,m=3,
    title="System matrix with boundary conditions") 

b[idx] .= 0;                 # Dirichlet values

u = A\b;
U = unvec(u)

• Accuracy and efficiency assuming square grid and finite differencing
  • \(\mathbf{A}\) has size \(N = O(n^2)\) with sparse Cholesky \(O(n^4)\)
  • \(n\) increases \(\rightarrow\) lower truncation error, more operations
  • for fixed run time, \(T\), convergence is \(O(1/\sqrt{T})\) where \(n = O(T^{1/4})\)
  • for dense Chebyshev spectral discretization, convergence is \(O(K^{-T^{1/6}})\) or \(O(K^{-n})\)