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})\)