Hello, I have been trying to solve vorticity and streamfunctions using finite-difference methods using linear algebra.
However, when I calculate the streamfunction \psi using the formula \zeta = \nabla^2\psi
I get extremely, extremely weird results.

function createΔmat!(Δmat::Array{<:Integer,2}, nx::Integer, ny::Integer)
ipt = 0; npts = nx * ny;
for iy = 1 : ny, ix = 1 : nx; ipt += 1
Δmat[ipt,ipt] = -4
ixb = ix - 1; if ixb == 0; ixb = nx end
ixa = ix + 1; if ixa > nx; ixa = 1 end
i1 = (iy-1) * nx + ixb
i2 = (iy-1) * nx + ixa
i3 = (iy-2) * nx + ix; if i3 < 0; i3 = mod(i3,npts); elseif i3 == 0; i3 = npts end
i4 = iy * nx + ix; if i4 > npts; i4 -= npts end
Δmat[i1,ipt] = 1
Δmat[i2,ipt] = 1
Δmat[i3,ipt] = 1
Δmat[i4,ipt] = 1
end
return
end
function findψfromq(q::Vector{<:Real},Δmat::AbstractArray;f::Real=0)
return Δmat \ (q .- f)
end

Looks like the matrix is singular. The smallest singular value is 0 - I get

julia> minimum(svd(Δmat).S)
4.601377645091019e-16

The matrix has a null space spanned by the vector of all ones, which you can see by

julia> norm(sum(Δmat,dims=2))
0.0

This is a pretty well-known feature of matrices resulting from discretizing the Laplacian with pure Neumann boundary conditions.

If you want to solve this system, you probably need to either use a pseudo-inverse or project the null space out of the solution. The former is expensive, while the latter would make Δmat dense unless you did this projection in each iteration within an iterative solver for matrix equations.

Edit: as @Paulo_Jabardo points out, you can also augment the system with constraints to remove the null space. This can be done using Lagrange multipliers, e.g.,

e = ones(size(Δmat,2))
b = e/sum(e)
A = [Δmat b;
transpose(b) 0]
b = [qvec
0]
u = A\b
ψ = u[1:end-1] # discard last entry = Lagrange multiplier enforcing a zero mean constraint

Yeah, I really should’ve mentioned that. Basically setting any constraint which specifies psi up to a constant should work (setting at a single node, setting the average, etc).

Thank you guys so much for your suggestions! I used the pinv() function, which as you mentioned is computationally expensive, but once I’m done with this problem set I’ll look into the methods you suggested to try and speed up the calculations.

It is very simple, setup yhe matrix, select the first row and make every off diagonal element zero and make the diagonal element 1. On right hand side, set the corresponding value to zero