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.
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.
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.
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.
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.
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).
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?
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):
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.
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".
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`.
Yes I know. It’s a solve operator…it’s equivalent to division only in the infinite precision world.
The \ operator differs from the / operator in that it doesn’t compute an inverse … it solves the system of equations. Solver algorithms are more numerically stable ( in that you’re much less likely to have large errors due to wacky input data).
I'd just say that solving the system of equations is the best way to divide by a matrix — that's why I put air quotes around "divide" above. In Julia, right-dividing matrices (with `A/B`) actually does the smart adjoint-commuting thing to left-divide it and do the solve (with `(B'\A')'`).
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...
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!
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.
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:
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!)
I agree with your point that it may be valuable to the reader to also mention common libraries or built-ins that solve the same problem that was discussed in the article somewhere near the end. More generally, linking to sources, docs, related work, or next steps can be very nice for readers to further their education. Framing your original comments with "the author should also mention... a standard solver like CHOLMOD" or "a standard solver like CHOLMOD ... as a sanity check that your implementation is working" as you described would have had a good chance of preventing our communication misfire and avoiding my (perhaps unnecessarily strong) snark.
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).
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.
