[enriquto]

Since you have written it as a symmetric, positive definite, sparse linear system, why don't use a standard solver like CHOLMOD which is available in Julia? (and behind octave's anti-slash operator). It should be faster than the ad-hoc single-scale Gauss-Seidel.

[yngvizzle]

While CHOLMOD is great, you cannot always use the Cholesky factorisation when you solve PDEs. For real-world simulations, we often have to solve systems with hundreds of millions, if not billions, of equations and in then case, even a highly optimised direct solver like CHOLMOD fails. The fill-in simply becomes too large.

For these small test cases, however, simply using CHOLMOD (or any other sparse solver) would do the trick perfectly.

[phkahler]

At the end the author mentions there are a large number of methods available for solving. Also mentions there are a much larger set of applications that those discussed.

[enriquto]

Sure! But it's a bit surprising that they do not use the language-provided linear solver and write simply f=A\b

[SonicScrub]

The goal of the article is educate on the mathematics, not a tutorial on how best to do mathematical modelling with Julia. An article like that is better served explaining the inner workings of the Black Box, rather than just using the Black Box.

[phkahler]

But they should mention the black box after the educational part.

[SonicScrub]

Why? The goal of the article is not to be a Julia tutorial. Julia is just a convenient demonstrator. Better to make the code as language independent as possible.

[icegreentea2]

Wait, how does this work? I'm really rusty, what would the f, A and b correspond to here?

[enriquto]

It's the common way to solve a linear system in octave, matlab and julia.

You have an invertible square matrix A, a vector b of the same dimension, and you want to find a vector x such that "A*x=b". Then you write "x=A\b", which is like "x=A^(-1)*b" but does not get to compute the full inverse matrix (which is useless).

[icegreentea2]

Sorry, yes I know about the syntax. I'm just struggling with what exactly you would be plugging in. Like with respect to any of the example problems given, what would A and x and B be?

[mdturnerphys]

You create matrices that represent the operator. For instance, suppose you had a 1D function and were interested in the evenly spaced points (x1, x2, x3, x4, x5), with the function values at these points being

  y1
  y2
  y3
  y4
  y5

If the differential equation had a first derivative in it, you'd construct something like this, multiplied by some constant (e.g. 1/(2*dx) for an unscaled derivative):

  ? ? ? ? ?
  -1 0 1 0 0
  0 -1 0 1 0
  0 0 -1 0 1
  ? ? ? ? ?

So the derivative at each element is defined by the difference of the next and previous elements. Multiplying the column of function values by this gives you the derivatives. This doesn't work for the first and last element, and in fact you'll usually modify these rows depending on what boundary condition is needed, so I've just left them filled in with "?".

For a solver, you don't know what the y values actually are, so you construct a column that corresponds to the right side of the differential equation. For instance, if the equation was something like the trivial dy/dx = c and you were using the operator above the column would be

  ?
  c
  c
  c
  ?

with the first and last values to be filled in based on the boundary conditions. You then left matrix divide that by the operator matrix (i.e. multiply it by the inverse of the operator matrix). That gives the solution to the equation.

This is just a simple example and in practice the matrices will be larger and more complex.

[enriquto]

x is f in the article

A is the five-point laplacian (a symmetric matrix with five non-zero diagonals), as implemented in the "step" function in the article, let's call it L

b is the datum h

If you set-up the sparse matrix L and the vector b (with the same code that performs the iterations in the article), then you can solve Poisson equation "L*f=h" by doing "f=L\h".

[hasmanean]

Don’t forget APL. I can’t say for sure but I imagine this \ operator came from there.

[mbauman]

It's not some obscure symbol; it's just division. It just so happens that we typically "divide" matrices from the left when solving equations like this (and matrix "division" isn't commutative) so instead of `a/b` it's `b\a`.

[ChrisRackauckas]

    using LinearAlgebra, SparseArrays
    N = 10
    D1 = sparse(Tridiagonal(ones(N-1),-2ones(N),ones(N-1)))
    Id = sparse(I,N,N)
    A = kron(D1,Id) + kron(Id,D1)

You just use the Kronecker product to build the matrix in the different directions. This, along with the relationship to convolutional neural networks, is described in the MIT 18.337 course notes: https://mitmath.github.io/18337/lecture14/pdes_and_convoluti...

[infogulch]

Wow! It seems you have the magical capacity to ingest the reference to an equation and instantly derive an intuition for how it works and what it's useful for. Learning technology this advanced has never been seen before by humankind, I hope you share it with the rest of us!

[burnished]

Everything being referenced is stuff you'd learn in a typical technical education as part of either calculus or statistics. Your reply comes off as of defensive in a way that implies that you're shocked some one would know any of this stuff.

[infogulch]

My reply stems from enriquto's misunderstanding of the purpose of the article, which is the "typical technical education" itself. It's like they are wondering why the article even exists, and isn't just a one line reference to the julia docs. Clearly there's nothing wrong with already having specific knowledge of a subject, but questioning the purpose of technical education because you already have it is bizarre.

Maybe an analogy would better explain my perspective. I imagine that enriquto would greatly appreciate my latest article, reproduced in its entirety below:

# Learn how to write a JSON parser

> j = JSON.parse("[1,2]")

Fin.

[enriquto]

Sorry for the misunderstanding, then. It was not at all my purpose to disparage this article. It is a lovely article and very clearly written and illustrated.

I'll state my point following your json parser example. If you write an article about the implementation of several json parsers, you may still want to call JSON.parse at the end, as a sanity check that your implementation is working. The function is right there and you may as well say that!

In the present case, since the author has already set-up explicitly Poisson equation as a linear system of equations, it would make sense to call julia's built-in solver. (If only to marvel that it is much, much faster than the simple methods shown before, thus it must make some really fancy stuff inside!)

[alisonkisk]

Please stop attacking other commenters. instead of attacking (the same person, twice!), try to understand what is being said, or ask questions.

[cpp_frog]

Numerical solutions of PDE are hardly studied in a typical technical education in calculus (it requires more theoretical machinery) or statistics (PDEs in statistics is somewhat narrow and specialized).

[Porygon]

Conjugate gradient descent with a multigrid preconditioner also works quite well in my experience, especially for larger systems.

[slavik81]

I've never worked with multi-grid, but I assume the preconditioner is also based on the Cholesky factorization? Incomplete Cholesky was pretty effective as the preconditioner for the pressure solve in my toy fluid sim.