Finite Difference, calculating the differentiation matrix, etc

Specific Domains Numerics
1

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

the construction of q (vorticity) is as follows:

dx   = 0.1
xvec = -1:dx:1; nx = length(xvec)
yvec = -1:dx:1; ny = length(yvec)
qvec = exp.(-(xvec.^2 .+ (yvec'/2).^2))[:]; npts = nx * ny
ψvec = zeros(npts)
Δmat  = zeros(Int8,npts,npts); createΔmat!(Δmat,nx,ny);

Does anyone know why the below result would occur? This really doesn’t make much sense to me

julia> findψfromq(qvec,Δmat)
 441-element Array{Float64,1}:
 -8.274473586102755e16
 -8.27447358610275e16
 -8.274473586102747e16
 -8.274473586102758e16
 -8.27447358610276e16
  ⋮
 -8.274473586102747e16
 -8.274473586102763e16
 -8.274473586102758e16
2

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

If I do this, I get the following plot from

import VectorizedRoutines.Matlab: meshgrid
x,y = meshgrid(xvec,yvec)
scatter(vec.((x,y,ψ)),zcolor=ψ)

Screen Shot 2020-09-24 at 7.30.28 PM
Screen Shot 2020-09-24 at 7.30.28 PM1062×716 222 KB

2 Likes
3

Wouldn’t setting the value of psi at a single node equal to zero solve the problem?

2 Likes
4

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).

5

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.

6

This method works very well and gives very similar results to using pinv!!! Thank you so much!!!

1 Like
7

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