# Finite Difference, calculating the differentiation matrix, etc

**URL:** <https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180>\
**Category:** Numerics\
**Created:** [September 24, 2020, 5:59am UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180 "2020-09-24T05:59:04Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![natgeo-wong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/natgeo-wong/32/10965_2.png) [@natgeo-wong](https://discourse.julialang.org/u/natgeo-wong)\
**Post date:** [September 24, 2020, 5:59am UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180/1 "2020-09-24T05:59:05Z")

</div>

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.

```julia
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:

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

```

---

<div class="post-metadata">

**Author:** ![jlchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlchan/32/10958_2.png) [@jlchan](https://discourse.julialang.org/u/jlchan)\
**Post date:** [September 24, 2020, 6:44pm UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180/2 "2020-09-24T18:44:33Z")

</div>

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

```julia
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
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.,

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

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

```

 ![Screen Shot 2020-09-24 at 7.30.28 PM](https://global.discourse-cdn.com/julialang/original/3X/1/9/1971d0012ee6cb4fd9cf42aba32821aca332e58b.png)

---

<div class="post-metadata">

**Author:** ![Paulo\_Jabardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paulo_jabardo/32/3196_2.png) [@Paulo\_Jabardo](https://discourse.julialang.org/u/Paulo_Jabardo)\
**Post date:** [September 25, 2020, 12:12am UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180/3 "2020-09-25T00:12:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jlchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlchan/32/10958_2.png) [@jlchan](https://discourse.julialang.org/u/jlchan)\
**Post date:** [September 25, 2020, 12:22am UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180/4 "2020-09-25T00:22:16Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![natgeo-wong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/natgeo-wong/32/10965_2.png) [@natgeo-wong](https://discourse.julialang.org/u/natgeo-wong)\
**Post date:** [September 25, 2020, 5:20am UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180/5 "2020-09-25T05:20:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![natgeo-wong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/natgeo-wong/32/10965_2.png) [@natgeo-wong](https://discourse.julialang.org/u/natgeo-wong)\
**Post date:** [September 25, 2020, 3:56pm UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180/6 "2020-09-25T15:56:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Paulo\_Jabardo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paulo_jabardo/32/3196_2.png) [@Paulo\_Jabardo](https://discourse.julialang.org/u/Paulo_Jabardo)\
**Post date:** [September 25, 2020, 4:01pm UTC](https://discourse.julialang.org/t/finite-difference-calculating-the-differentiation-matrix-etc/47180/7 "2020-09-25T16:01:43Z")

</div>

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