7 ms·

Poisson's Equation

766 points by mferraro89 5y ago

heinrichhartman 5y ago
This guy has exactly two posts on his blog. The other one is on a completely different topic but great read as well: https://mattferraro.dev/posts/cnc-router https://mattferraro.dev/posts/cnc-router
enriquto 5y ago
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.
- phkahler 5y ago
  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 5y ago
    Sure! But it's a bit surprising that they do not use the language-provided linear solver and write simply f=A\b
    - infogulch 5y ago
      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 5y ago
        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.
        cpp_frog 5y ago
        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).
        infogulch 5y ago
        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 5y ago
        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!)
        5y ago
    - SonicScrub 5y ago
      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 5y ago
        But they should mention the black box after the educational part.
        SonicScrub 5y ago
        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 5y ago
      Wait, how does this work? I'm really rusty, what would the f, A and b correspond to here?
      - enriquto 5y ago
        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).
        hasmanean 5y ago
        Don’t forget APL. I can’t say for sure but I imagine this \ operator came from there.
        mbauman 5y ago
        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`.
        hasmanean 5y ago
        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).
        mbauman 5y ago
        Yes, I know. :) 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')'`).
    - deleted 5y ago
      [deleted]
- Porygon 5y ago
  Conjugate gradient descent with a multigrid preconditioner also works quite well in my experience, especially for larger systems.
  - slavik81 5y ago
    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.
- yngvizzle 5y ago
  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.
zitterbewegung 5y ago
Sorry but can there be more context to why it is a powerful tool?
- wildmanx 5y ago
  This is a typical issue with HN posts. Some poor soul wrote a somewhat competent and maybe even lengthy blog post / article about something they really care about and are knowledgeable about. It may be directed at a specific audience, or maybe just screaming into the void to record down some insight they had for themselves to read again later, or similar. And they use a more-than-necessary general title like "the best tool you'll ever see" with "you" meaning either just themselves or a narrow target audience or so. And then somebody comes along who finds it interesting, submits it to HN, it makes front page, and now it looks like the poor author with their more-general-than-needed title is making a general statement about sth for the changed audience which is the HN crowd, but which is distinctly different from the original target audience. Case in point: The articles author self describes as: "I'm an aerospace engineer that writes software. I love math and science, and I have two cats." For an aerospace engineer this all makes a lot of sense and is a super great tool, I'm sure. It's just not for the overall HN crowd. And it's not the authors fault.
  - rawtxapp 5y ago
    Just to point out, it seems like the author themselves submitted the story to HN. Still I agree with your reasoning, but I think slight clickbaity titles do get clicks which is why we keep seeing them.
    - MontyCarloHall 5y ago
      Yup. This wouldn’t get nearly the same amount of attention if the title were “Poisson’s equation is the most powerful tool in your toolbox for finding steady-state solutions to the heat equation with arbitrarily placed sources”
      - crazygringo 5y ago
pgustafs 5y ago
Great post, one nitpick -- I wouldn't say that a matrix is a "sparsely defined" function, but rather a function defined on a finite grid. It might also be worth pointing out that same approach works for any graph, not just a grid.
- wildmanx 5y ago
  Also, what's confusing is that algebra usually uses matrices to describe linear functions from n-dimensional to m-dimensional vector spaces. Matrix has n rows, m columns, you give it an n-dim vector and after matrix multiplication you get back an m-dim vector. The author uses a matrix quite differently. You give it two integer coordinates i and j and it gives you the value at position (i, j) back. That's a valid use, but not quite what you'd expect in a math-oriented article.
  - burnished 5y ago
    Can you link to context for this? I learned both in linear algebra, so it seems like either would be just as 'expected'.
    - mixedmath 5y ago
      Here's a concrete example. The first matrix in the post is f = [[1, 1, 1], [1, 1, 1], [1, 1, 1]]. In linear algebra, we would interpret this as a linear map. A true equation would be f([1, 2, 3]^T) = [6, 6, 6]^T (where I'm using ^T to mean "transpose to a column vector"). But here, the author means f(1, 2) = 1, i.e. the (1,2) coordinate of the matrix is 1.
      - burnished 5y ago
        Thank you! Yes, I agree, thank you for explaining that to me.
        wildmanx 5y ago
        And interestingly, both are connected. If d_i somewhat hand-wavingly expresses the vector d_i = (0, ..., 0, 1, 0, ... 0) with the 1 at position i, then given matrix M you can do f(i, j) := d_i^T * M * d_j The RHS is using classical matrix multiplication, and the function value will be the matrix' entry at column i, row j.
pm90 5y ago
This was a very interesting read even as someone who probably has no practical use for these tools.
- davidkuhta 5y ago
  > probably can lead you to some fun places.
  - xbar 5y ago
    I am intrigued, but I think I still need a bridge to help me get to practical uses in my my fields involving (largely) non-physical systems.
jbay808 5y ago
A friend of mine broke a badminton racket during a match, and I was struck by how the sharply bent and twisted metal rim was transformed into a smooth, continuously double-curved surface by the racket weave. I looked closely at the balance of tension in the woven cord, thought of how it resembles Poisson's equation, and suddenly it all made sense. Edit - it looked something like this: https://thumbs.dreamstime.com/b/broken-badminton-racket-photo-taken-malaysia-45327300.jpg https://thumbs.dreamstime.com/b/broken-badminton-racket-phot...
- codethief 5y ago
  Hmmm, this looks more like a minimal surface, i.e. a solution to the minimal-surface equation[0], than a solution to Poisson's equation. Then again, both equations are of elliptic type. Some links for people who've never heard of minimal surfaces: https://en.wikipedia.org/wiki/Minimal_surface https://en.wikipedia.org/wiki/Minimal_surface https://minimalsurfaces.blog/ https://minimalsurfaces.blog/ (lots of illustrations) https://makmanx.github.io/math3435s18/talks/MSE.pdf https://makmanx.github.io/math3435s18/talks/MSE.pdf (brief intro with historical remarks and illustrations) [0]: More specifically, it's a solution to Plateau's problem: https://en.wikipedia.org/wiki/Plateau%27s_problem https://en.wikipedia.org/wiki/Plateau%27s_problem
  - jbay808 5y ago
    Isn't Poisson's equation basically describing a minimal surface for small z? I'm not saying the badminton racket follows exactly a (discrete) 2D Poisson equation. But it's certainly related enough to be more than a surface similarly. The cords are under high tension, which means that any curvature along x (that is, dz^2/dx^2) will result in a net z-axis tension force unless balanced by an oppositely curved cord running in the y direction. Since it's in static equilibrium, there can be no unbalanced forces and so that must be the case. Therefore at each intersection, (d^2/dx^2 + d^2/dy^2)z = 0, which is Poisson's equation in 2D for z height being the function. Approximately, assuming equal tension in x and y, small z, and so on.
    - codethief 5y ago
Jeff_Brown 5y ago
EDIT: I'm leaving this here to help anyone else who might have been confused by this, which I imagine is likely. What confused me is that the author is not treating the matrix as a function from vectors to vectors, as is the customary way to treat matrices as functions. Rather, they're using the matrix to represent a sparse, regular sampling of a function from vectors to scalars. --- This article makes no sense right off the bat. Here's the first substantive passage: "[Laplace's Equation means] Find me a function f where every value everywhere is the average of the values around it ... In this post, when we talk about a function f we mean a 2D matrix where each element is some scalar value like temperature or pressure or electric potential ... If it seems weird to call a matrix a function, just remember that all matrices map input coordinates (i,j) to output values f(i,j). Matrices are functions that are just sparsely defined. This particular matrix does satisfy Laplace's equation because each element [of the matrix] is equal to the average of its neighbors." The values of a function are the outputs it maps its inputs to. The elements of the matrix are neither inputs nor outputs.
- freeone3000 5y ago
  The matrix is the function. The elements are the outputs. The coordinates are the inputs. Take, for example, f(x) = x*x. Its matrix would be: f = [0, 1, 4, 9, 16].
- techas 5y ago
  A matrix can be seen as the discrete representation of a function…
- gmmeyer 5y ago
  Poisson's Equation takes another function as an input, the matrix is the representation of the output of said function
- jungturk 5y ago
  Think of the matrix as a precomputed lookup table. Given the arguments to the function, locate the cell in the matrix and use its value as the result of the function.
prionassembly 5y ago
My (unorthodox and somewhat rickety) note-taking gizmo uses Poisson's equation to classify (continuously) entries. Basically the note-taking gizmo is a graph. Nodes are given conceptual masses either through pagerank or betweenness centrality (i.e. either through how many random walks or how many shortest paths cross a node). Then we calculate a potential energy (gravity potential) if we by inverting the graph laplacian (a few methods are available). Special attention is given to nodes that "float the most. E: forgot to link to it! https://github.com/asemic-horizon/sursis/ https://github.com/asemic-horizon/sursis/
- chris_st 5y ago
  That looks really cool -- I presume it has some way to enter more than a single word/phrase in a node? Not sure linking individual words is useful for me :-) I should just try it, of course...
- throwamon 5y ago
  I guess I don't have the background to see how useful this can be, but something tells me it can be very useful. Would you mind doing an ELI5?
mixedmath 5y ago
I think this is a beautiful article. There's code, there's math. There are many plotted examples using a variety of different plotting techniques (in total, representing a 2d array as data, or in black/white, or with a terrible 'jet' colormap, or as a 3d terrain. I very much appreciate this sort of post.
- nickponline 5y ago
  Maybe you'll like these: https://nickp.svbtle.com/ https://nickp.svbtle.com/ (shameless plug)
  - lesquivemeau 5y ago
    Thanks for this
pphysch 5y ago
The notion that the simple Laplace solver can scale to grids of arbitrary size, without modification, is a bit misleading for practical purposes. Computational performance will tank or zero out if memory hierarchy constraints are not considered. The author does mention high-performance solutions like multigrid. However, even a basic successive overrelaxation algorithm like the one shown can be partitioned and parallelized, and it is a very good programming exercise to implement a partitioning scheme using MPI or even a low-performance messaging library (or even just optimize for cache sizes on a single device, with no network transport). Like the transition from the elegant, 5-character Laplace equation to the relatively verbose and complex numerical Julia solver, there is an additional and necessary step in making the numerical solution further scalable with present technology. In particular, the notion that the computational boundaries map nicely to the physical boundaries must be thrown out, because now we must respect the layer of "virtual boundaries" between the partitions.
thendrill 5y ago
Very very well written and enlightening article. Love this guy's writing.
seemslegit 5y ago
That's such a presumptuous title, the author does not know what other powerful tools my toolbox lacks.
SavantIdiot 5y ago
I could have used this 32(!) years ago when I was struggling in college. (This and 3b1b.) It amazes me just how many key topics were so inaccessible to the majority of the class at engineering school. I base this on observations from group study sessions and the hyper-aggressive test curves. I knew lots of people who never got the hang of div/grad/curl, or a Jacobians, or eignenvectors, or Z-transforms... These are key engineering concepts, you'd think colleges would bend over backwards to make sure these concepts are learned as succinctly as possible rather than add a curve to a test that makes a 23 out of 100 an "A" grade. I'm digressing, and complaining, but the counter argument has always been: you're not supposed to learn everything in college, you're supposed to learn how to learn. Sure, right, but who has time to keep learning advanced calculus after college? (Well, I still study math & physics for fun, but over the course of decades, not years.) Not being able to see the world through these lenses I think means missing key engineering perspectives and relationships. Anyway, very well written article.
- lordnacho 5y ago
  I feel the same. How can it be that I went to a world famous institution providing 2-to-1 student-teacher ratios, but I still think the best explanations are these modern internet explanations? I guess the best explanations just bubble up in the modern environment. > you're not supposed to learn everything in college, you're supposed to learn how to learn But to learn how to learn, you gotta learn some things to a somewhat decent degree. I think at some point you need to have these linalg/divgradcurl things down, if only briefly. You might forget any particular topic, but if you've indexed it you should be able to pick it up again, particularly in the modern learning environment. Just imagine coding without access to StackOverflow.
  - pfortuny 5y ago
    Honest question. 2-to-1 student-teacher ratio?? Where is that? Totally stunned.
    - ectopod 5y ago
      Oxford tutorials do this. They typically have one, two or three students.
panic_on_oops 5y ago
Wonderful post, thank you OP!
tobmlt 5y ago
See also from Keenan Crane and company (discrete differential geometry): http://ddg.cs.columbia.edu/SGP2014/LaplaceBeltrami.pdf http://ddg.cs.columbia.edu/SGP2014/LaplaceBeltrami.pdf “The Swiss Army knife of geometric operators.” I always thought that was cool since I usually think of diffusion in the context of fluid flow.
niffydroid 5y ago
You lost me at equation
Chris2048 5y ago
Reminds me of 3Blue1Brown "Divergence and curl": https://www.youtube.com/watch?v=rB83DpBJQsE https://www.youtube.com/watch?v=rB83DpBJQsE
s-macke 5y ago
You can also use it to solve labyrinths. Just put a high pressure at the beginning and a low pressure at the end. Solve the Poisson equation. The path through the labyrinth is always the steepest slope. In [1] you can see a small implementation of the idea. [1] https://simulationcorner.net/maze/ https://simulationcorner.net/maze/
- setr 5y ago
  More generally, I believe this is exactly the same as flow-field pathfinding, which conveniently has the property that you end up with a solution to the whole map -- so you've solved all pathfinding for all entities anywhere in the map, in one shot. Great for hosting a ridiculous amount of entities on a static (or slow-changing) map. Brogue's creator also "invented" djikstra maps, which I believe is also exactly the same[0][1], to handle AI strategic pathfinding (e.g. avoid hazards while reaching treasure). In Game AI Pro (I forget which article/book), someone suggests taking AI preferences (e.g. hunger, health) and computing a flow-field per preference (e.g. a food-map, a danger-map, etc), and multiplying the values against the preference to act as a weight (so food-map * %hunger, danger-map * %health), and summing the maps together to produce the final map used for pathfinding. Notably, the food-map can be shared by all entities consuming the same kind of food -- only the weight has to be re-calculated, and the final sum. [0] http://www.roguebasin.com/index.php?title=The_Incredible_Power_of_Dijkstra_Maps http://www.roguebasin.com/index.php?title=The_Incredible_Pow... [1] http://www.roguebasin.com/index.php/Dijkstra_Maps_Visualized http://www.roguebasin.com/index.php/Dijkstra_Maps_Visualized
alisonkisk 5y ago
The article never explains why ∇2 means "average of neighbors". It's the divergence of the gradient, which is (one kind of) n-dimensional 2nd derivative. In a one-dimensional function, the second derivative is 0 when there is no curvature, aka a straight line (of any slope), and any point on a line is equal to the average of its neighborhood. A plane also has this property, but in 2+ dimensions you can also make other shapes (like saddles), that are made up of lines (like a plane) but the lines are twisted relative to each other in interesting ways (like "string art"). You can also visualize (aka impose a coordinate system for) these surfaces as having positive curvature (concave up) in one direction, and exactly opposite negative curvature (convex up, or concave down) in the orthogonal direction, summing to 0.
vlmutolo 5y ago
> It is customary when simulating heat flow to use a wacky color palette where red is hot and blue is cold, with all kinds of intermediate colors in between In an otherwise excellent article, this is the only issue I could find. We really need to stop using/recommending/normalizing rainbow color maps (i.e. jet). They actively confuse readers by creating visual artifacts that aren't actually in the data. This article has some great explanations and visuals. https://jakevdp.github.io/blog/2014/10/16/how-bad-is-your-colormap/ https://jakevdp.github.io/blog/2014/10/16/how-bad-is-your-co... The original post uses a rainbow color map to represent temperature-related things because having a diverging color map is often a useful intuition for temperature heat maps. But in that case, we should prefer one of the following diverging color maps listed on the matplotlib site (this list definitely isn't exhaustive, but it is helpful). https://matplotlib.org/stable/_images/sphx_glr_colormaps_004.png https://matplotlib.org/stable/_images/sphx_glr_colormaps_004... More on matplotlib's well-chosen color maps: https://matplotlib.org/stable/tutorials/colors/colormaps.html https://matplotlib.org/stable/tutorials/colors/colormaps.htm...
- just_temp 5y ago
  Can not agree with this more, if people want to plot something that is linear please use a perceptually linear colormap! Just a one second glance at the Mona Lisa in rainbow/Jet is enough to make you gouge your eyes out. https://peterjamesthomas.com/2017/09/15/hurricanes-and-data-visualisation-part-ib-the-mona-lisa/ https://peterjamesthomas.com/2017/09/15/hurricanes-and-data-... For a more technical description the information behind the newer matplotlib defaults, particularly the scipy talk, is great. https://bids.github.io/colormap/ https://bids.github.io/colormap/ And for those that do not like the matplotlib options, colorcet provides a wider range of alternatives that are not trash (unlike Jet) https://colorcet.holoviz.org/index.html https://colorcet.holoviz.org/index.html
  - vlmutolo 5y ago
    Thanks for the links! I especially like the Mona Lisa example. I'll probably steal that the next time this topic comes up.
cornel_io 5y ago
> This oval (called a separatrix) has the special property that no wind flows through it at all. It acts just like a solid surface. We have already assembled a reasonably accurate simulation of how air flows around an ellipse in some confined space like a wind tunnel! My understanding of airflow simulation (undergrad-level at best) is that the correct boundary condition is almost without exception the no-slip one: air should be stationary at the surface of each object, not just have zero flow across it. Am I correct that the calculation mentioned above really only applies to the "dry liquid" scenario where there is no drag and zero viscosity?
- mferraro89 5y ago
  You are correct! The simulation as written does not handle the no slip condition. The simulated "air" is inviscid and irrotational. You would need to tackle a few more things in order to have a really accurate simulation.
  - SubiculumCode 5y ago
    Thanks for the awesome write-up.writer's note. The first paragraph of your conclusion would have served better as your introductory paragraph.
sgarrity 5y ago
I thought poisson's equation was expressed as: <><
toolslive 5y ago
I consider Pagerank to be a discrete variation of the poission equation.
Synaesthesia 5y ago
We referred to this method of numerically solving Poisson's equation by successively averaging values, as the method of relaxation.
Lichtso 5y ago
Last year a Monte Carlo approach to estimate the solution quickly was discovered. It works somewhat similar to diffusion curves. https://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/index.html https://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryP...
WalterBright 5y ago
A very understandable explanation of it, and its uses.
Robotbeat 5y ago
I have used Poisson’s Equation recently. Actually, something similar. Been coding up some thermal conductivity calculations form scratch. What’s nice about the Laplace equation is that there are some exact analytical solutions in some situations which is really helpful for validating numerical codes and sanity checking. It’s also simple enough to implement in an Excel spreadsheet, Color coding the cells to indicate temperature. Good sort of “no code” example of the concept.
daleroberts 5y ago
Here is my implementation in python of the Poisson equation on an arbitrary 2D domain using the finite element method. I used this for teaching a course in partial differential equations: https://github.com/daleroberts/poisson https://github.com/daleroberts/poisson
MauranKilom 5y ago
> From here we could use Bernoulli's equation to find the pressure distribution, which we could integrate over the surface to find drag and lift and so on. With a few tweaks we could simulate rotational flow, vortex panels, real wing profiles, and so on. > With just a few simple building block we're already edging up on real computational fluid dynamics. All this just by adding up some matrices! Correct me if I'm wrong, but the only "real" CFD you could solve with this are incompressible potential flows [0]. Solving Navier Stokes is clearly not just "a few tweaks away" from the Laplace equation, but I would be curious which tweaks would take you to e.g. rotational flows. [0]: https://en.wikipedia.org/wiki/Potential_flow https://en.wikipedia.org/wiki/Potential_flow
powderpig 5y ago
I too hate the coloured plots but I can tell you as a Stress & Structures engineer who works with ANSYS regularly, it does help when assessing material limits and strain energies. This is especially true when you're working with models where temperature is time dependent. A great article, one to bookmark for sure.
woopwoop 5y ago
Why is the laplacian so ubiquitous? Well, locally any reasonable PDE is well approximated by a linear one. And linear PDEs are of the form Lu = f for a linear differentiable operator L. But why does the particular case of L = laplacian show up so often in physics? Galilean invariance says that the laws of physics should be invariant under translation and rotation. Checking this against Lu = f, you can verify that one requires that L commute with translations and rotations. That is L(u_h) = (Lu)_h, where u_h(x) = u(x+h), and similarly for rotations. What you can show (pretty easily on the Fourier side) is that every linear differential operator with these properties is a polynomial in the Laplacian.
- amirkdv 5y ago
  > you can show [...] every linear differential operator [that commutes with translations and rotations] is a polynomial in the Laplacian. I may be tripping over something here but this doesn't sound right. If you mean polynomial in the real number sense, i.e. Lu = a0 + a1 Delta(u) + a2 Delta(u)^2 + ... (where Delta = Laplacian, and a0, a1, ... real numbers), then is this true? The famous wave operator doesn't have this form. And if you mean "polynomial" as in a series over function space, i.e. Lu ~ a0 + a1 Du + a2 D^2u + ... (infinite terms, not equality but convergence) where D is the usual differential operator and a0 is a number, a1 is a 1-d vector, a2 is a 2d matrix, so on, then this is standard calculus of variations on any reasonable function space. Taylor series for function spaces if you want. I don't think that's limited to nice Galilean operators.
  - woopwoop 5y ago
    Thanks for pointing this out, I should have been more precise. I was referring to time invariant equations. You are right that this doesn't apply to the wave operator (which also does not commute with rotations in R^4). But in space only, yes polynomial means polynomial with constant coefficients. The proof basically goes that commuting with translations implies immediately that L has constant coefficients. Then on the Fourier side applying L translates to multiplying by a polynomial, and commuting with rotations translates to the claim that that polynomial is rotation invariant. And every rotation invariant polynomial in several variables is p(|x|^2) for some polynomial p in on variable. Then on the spacial side p(|x|^2) translates to L = p(laplacian).
Attained 5y ago
What does it mean? I fully expected to explain the semantics and yet I still don't know what that triangle squared means. Oh well.
steve76 5y ago
This article details matrices as fixed known sets of data. Much more interesting is the vector calculus. Vectors as the result of the sum of their components are good. But think of the reverse. Your model outputs the current vector. Your model predicts the next vector in sequence. Just from current and next, your model provides any amount of underlying components, as high resolution as you want. Nature fails before your model does, through the lack of continuity in steps and impulses and the mean value theorem no longer holds. Vector analysis was pioneered by Josiah Gibbs, someone who doesn't nearly get enough credit. Up there with Boltzmann, Maxwell, Planck, and Einstein.
carodgers 5y ago
Beautifully done. Everything was clear and intelligible. I have a CS background, but no physics/aerodynamics background. What are the minimum additional steps beyond this tutorial which could produce an aerodynamically correct 2D wing simulator? (With turbulence, not a steady-state solution.) I've come back to tackle this topic intermittently but have never cracked it. To anyone with expertise here who could share an overview and links, I'd be grateful.
- btrettel 5y ago
  There are many different ways to do what you'd like. The easiest starting point would probably be this tutorial: https://github.com/barbagroup/CFDPython https://github.com/barbagroup/CFDPython But that won't handle turbulence. The real "turbulence problem" is that computing actual turbulent flows requires enormous computational resources. So instead of solving the Navier-Stokes equations, related equations with lower computational cost are solved. Because of how these equations are developed, they require modeling of "unclosed" terms, and this is a likely source of inaccuracy. If you want something relatively simple, you could take the RANS approach and use the Spalart-Allmaras model: https://www.cfd-online.com/Wiki/Introduction_to_turbulence/Reynolds_averaged_equations https://www.cfd-online.com/Wiki/Introduction_to_turbulence/R... https://www.cfd-online.com/Wiki/Spalart-Allmaras_model https://www.cfd-online.com/Wiki/Spalart-Allmaras_model How to implement the changes to the final part of Lorena Barba's tutorial should be fairly obvious by the time you get there.
  - carodgers 5y ago
    The CFD Online site is a great resource that I had missed before. Thank you for the links.
ybarhighbar 5y ago
Always looking to interrelate my knowledge that sometimes appears to more disparate that I care to consider. The simplicity of considering a matrix function in this way is very elegant. This type of methodology can help demystify some of the earlier steps in modeling by adding a simple intuitive picture that encourages a connection to calculation from the get go.
aabbcc1241 5y ago
The article explain clearly. Made a toy canvas to play with poisson https://github.com/beenotung/poisson-canvas https://github.com/beenotung/poisson-canvas
amai 5y ago
"you should know that Laplace's Equation has a close cousin called the Biharmonic Equation... Which is widely used in modeling quantum mechanics. " Can someone elaborate on that? The Wikipedia article (https://en.wikipedia.org/wiki/Biharmonic_equation https://en.wikipedia.org/wiki/Biharmonic_equation) does not mention anything about quantum mechanics use cases for the biharmonic equation.
vayhay 5y ago
Hiện tại hệ thống VayHay đang tuyển nhân viên Telesale khu vực cần thơ và các tỉnh thành khác toàn quốc làm việc 8h/ngày và có hỗ trợ vay tiêu dùng toàn quốc không thế chấp kì hạn và các gói vay linh hoạt, truy cập VayHay.vn để biết thêm thông tin

766 points by mferraro89 5y ago

This guy has exactly two posts on his blog. The other one is on a completely different topic but great read as well: https://mattferraro.dev/posts/cnc-router https://mattferraro.dev/posts/cnc-router

enriquto 5y ago

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.

phkahler 5y ago
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 5y ago
  Sure! But it's a bit surprising that they do not use the language-provided linear solver and write simply f=A\b
  - infogulch 5y ago
    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 5y ago
      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.
      - cpp_frog 5y ago
        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).
      - infogulch 5y ago
        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 5y ago
        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!)
        5y ago
  - SonicScrub 5y ago
    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 5y ago
      But they should mention the black box after the educational part.
      - SonicScrub 5y ago
        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 5y ago
    Wait, how does this work? I'm really rusty, what would the f, A and b correspond to here?
    - enriquto 5y ago
      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).
      - hasmanean 5y ago
        Don’t forget APL. I can’t say for sure but I imagine this \ operator came from there.
        mbauman 5y ago
        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`.
        hasmanean 5y ago
        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).
        mbauman 5y ago
        Yes, I know. :) 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')'`).
  - deleted 5y ago
    [deleted]
Porygon 5y ago
Conjugate gradient descent with a multigrid preconditioner also works quite well in my experience, especially for larger systems.
- slavik81 5y ago
  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.
yngvizzle 5y ago
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.

zitterbewegung 5y ago

Sorry but can there be more context to why it is a powerful tool?

wildmanx 5y ago
This is a typical issue with HN posts. Some poor soul wrote a somewhat competent and maybe even lengthy blog post / article about something they really care about and are knowledgeable about. It may be directed at a specific audience, or maybe just screaming into the void to record down some insight they had for themselves to read again later, or similar. And they use a more-than-necessary general title like "the best tool you'll ever see" with "you" meaning either just themselves or a narrow target audience or so. And then somebody comes along who finds it interesting, submits it to HN, it makes front page, and now it looks like the poor author with their more-general-than-needed title is making a general statement about sth for the changed audience which is the HN crowd, but which is distinctly different from the original target audience. Case in point: The articles author self describes as: "I'm an aerospace engineer that writes software. I love math and science, and I have two cats." For an aerospace engineer this all makes a lot of sense and is a super great tool, I'm sure. It's just not for the overall HN crowd. And it's not the authors fault.
- rawtxapp 5y ago
  Just to point out, it seems like the author themselves submitted the story to HN. Still I agree with your reasoning, but I think slight clickbaity titles do get clicks which is why we keep seeing them.
  - MontyCarloHall 5y ago
    Yup. This wouldn’t get nearly the same amount of attention if the title were “Poisson’s equation is the most powerful tool in your toolbox for finding steady-state solutions to the heat equation with arbitrarily placed sources”
    - crazygringo 5y ago

pgustafs 5y ago

Great post, one nitpick -- I wouldn't say that a matrix is a "sparsely defined" function, but rather a function defined on a finite grid. It might also be worth pointing out that same approach works for any graph, not just a grid.

wildmanx 5y ago
Also, what's confusing is that algebra usually uses matrices to describe linear functions from n-dimensional to m-dimensional vector spaces. Matrix has n rows, m columns, you give it an n-dim vector and after matrix multiplication you get back an m-dim vector. The author uses a matrix quite differently. You give it two integer coordinates i and j and it gives you the value at position (i, j) back. That's a valid use, but not quite what you'd expect in a math-oriented article.
- burnished 5y ago
  Can you link to context for this? I learned both in linear algebra, so it seems like either would be just as 'expected'.
  - mixedmath 5y ago
    Here's a concrete example. The first matrix in the post is f = [[1, 1, 1], [1, 1, 1], [1, 1, 1]]. In linear algebra, we would interpret this as a linear map. A true equation would be f([1, 2, 3]^T) = [6, 6, 6]^T (where I'm using ^T to mean "transpose to a column vector"). But here, the author means f(1, 2) = 1, i.e. the (1,2) coordinate of the matrix is 1.
    - burnished 5y ago
      Thank you! Yes, I agree, thank you for explaining that to me.
      - wildmanx 5y ago
        And interestingly, both are connected. If d_i somewhat hand-wavingly expresses the vector d_i = (0, ..., 0, 1, 0, ... 0) with the 1 at position i, then given matrix M you can do f(i, j) := d_i^T * M * d_j The RHS is using classical matrix multiplication, and the function value will be the matrix' entry at column i, row j.

pm90 5y ago

This was a very interesting read even as someone who probably has no practical use for these tools.

davidkuhta 5y ago
> probably can lead you to some fun places.
- xbar 5y ago
  I am intrigued, but I think I still need a bridge to help me get to practical uses in my my fields involving (largely) non-physical systems.

jbay808 5y ago

A friend of mine broke a badminton racket during a match, and I was struck by how the sharply bent and twisted metal rim was transformed into a smooth, continuously double-curved surface by the racket weave. I looked closely at the balance of tension in the woven cord, thought of how it resembles Poisson's equation, and suddenly it all made sense. Edit - it looked something like this: https://thumbs.dreamstime.com/b/broken-badminton-racket-photo-taken-malaysia-45327300.jpg https://thumbs.dreamstime.com/b/broken-badminton-racket-phot...

codethief 5y ago
Hmmm, this looks more like a minimal surface, i.e. a solution to the minimal-surface equation[0], than a solution to Poisson's equation. Then again, both equations are of elliptic type. Some links for people who've never heard of minimal surfaces: https://en.wikipedia.org/wiki/Minimal_surface https://en.wikipedia.org/wiki/Minimal_surface https://minimalsurfaces.blog/ https://minimalsurfaces.blog/ (lots of illustrations) https://makmanx.github.io/math3435s18/talks/MSE.pdf https://makmanx.github.io/math3435s18/talks/MSE.pdf (brief intro with historical remarks and illustrations) [0]: More specifically, it's a solution to Plateau's problem: https://en.wikipedia.org/wiki/Plateau%27s_problem https://en.wikipedia.org/wiki/Plateau%27s_problem
- jbay808 5y ago
  Isn't Poisson's equation basically describing a minimal surface for small z? I'm not saying the badminton racket follows exactly a (discrete) 2D Poisson equation. But it's certainly related enough to be more than a surface similarly. The cords are under high tension, which means that any curvature along x (that is, dz^2/dx^2) will result in a net z-axis tension force unless balanced by an oppositely curved cord running in the y direction. Since it's in static equilibrium, there can be no unbalanced forces and so that must be the case. Therefore at each intersection, (d^2/dx^2 + d^2/dy^2)z = 0, which is Poisson's equation in 2D for z height being the function. Approximately, assuming equal tension in x and y, small z, and so on.
  - codethief 5y ago

Jeff_Brown 5y ago

EDIT: I'm leaving this here to help anyone else who might have been confused by this, which I imagine is likely. What confused me is that the author is not treating the matrix as a function from vectors to vectors, as is the customary way to treat matrices as functions. Rather, they're using the matrix to represent a sparse, regular sampling of a function from vectors to scalars. --- This article makes no sense right off the bat. Here's the first substantive passage: "[Laplace's Equation means] Find me a function f where every value everywhere is the average of the values around it ... In this post, when we talk about a function f we mean a 2D matrix where each element is some scalar value like temperature or pressure or electric potential ... If it seems weird to call a matrix a function, just remember that all matrices map input coordinates (i,j) to output values f(i,j). Matrices are functions that are just sparsely defined. This particular matrix does satisfy Laplace's equation because each element [of the matrix] is equal to the average of its neighbors." The values of a function are the outputs it maps its inputs to. The elements of the matrix are neither inputs nor outputs.

freeone3000 5y ago
The matrix is the function. The elements are the outputs. The coordinates are the inputs. Take, for example, f(x) = x*x. Its matrix would be: f = [0, 1, 4, 9, 16].
techas 5y ago
A matrix can be seen as the discrete representation of a function…
gmmeyer 5y ago
Poisson's Equation takes another function as an input, the matrix is the representation of the output of said function
jungturk 5y ago
Think of the matrix as a precomputed lookup table. Given the arguments to the function, locate the cell in the matrix and use its value as the result of the function.

prionassembly 5y ago

My (unorthodox and somewhat rickety) note-taking gizmo uses Poisson's equation to classify (continuously) entries. Basically the note-taking gizmo is a graph. Nodes are given conceptual masses either through pagerank or betweenness centrality (i.e. either through how many random walks or how many shortest paths cross a node). Then we calculate a potential energy (gravity potential) if we by inverting the graph laplacian (a few methods are available). Special attention is given to nodes that "float the most. E: forgot to link to it! https://github.com/asemic-horizon/sursis/ https://github.com/asemic-horizon/sursis/

chris_st 5y ago
That looks really cool -- I presume it has some way to enter more than a single word/phrase in a node? Not sure linking individual words is useful for me :-) I should just try it, of course...
throwamon 5y ago
I guess I don't have the background to see how useful this can be, but something tells me it can be very useful. Would you mind doing an ELI5?

mixedmath 5y ago

I think this is a beautiful article. There's code, there's math. There are many plotted examples using a variety of different plotting techniques (in total, representing a 2d array as data, or in black/white, or with a terrible 'jet' colormap, or as a 3d terrain. I very much appreciate this sort of post.

nickponline 5y ago
Maybe you'll like these: https://nickp.svbtle.com/ https://nickp.svbtle.com/ (shameless plug)
- lesquivemeau 5y ago
  Thanks for this

pphysch 5y ago

The notion that the simple Laplace solver can scale to grids of arbitrary size, without modification, is a bit misleading for practical purposes. Computational performance will tank or zero out if memory hierarchy constraints are not considered. The author does mention high-performance solutions like multigrid. However, even a basic successive overrelaxation algorithm like the one shown can be partitioned and parallelized, and it is a very good programming exercise to implement a partitioning scheme using MPI or even a low-performance messaging library (or even just optimize for cache sizes on a single device, with no network transport). Like the transition from the elegant, 5-character Laplace equation to the relatively verbose and complex numerical Julia solver, there is an additional and necessary step in making the numerical solution further scalable with present technology. In particular, the notion that the computational boundaries map nicely to the physical boundaries must be thrown out, because now we must respect the layer of "virtual boundaries" between the partitions.

thendrill 5y ago

Very very well written and enlightening article. Love this guy's writing.

seemslegit 5y ago

That's such a presumptuous title, the author does not know what other powerful tools my toolbox lacks.

SavantIdiot 5y ago

I could have used this 32(!) years ago when I was struggling in college. (This and 3b1b.) It amazes me just how many key topics were so inaccessible to the majority of the class at engineering school. I base this on observations from group study sessions and the hyper-aggressive test curves. I knew lots of people who never got the hang of div/grad/curl, or a Jacobians, or eignenvectors, or Z-transforms... These are key engineering concepts, you'd think colleges would bend over backwards to make sure these concepts are learned as succinctly as possible rather than add a curve to a test that makes a 23 out of 100 an "A" grade. I'm digressing, and complaining, but the counter argument has always been: you're not supposed to learn everything in college, you're supposed to learn how to learn. Sure, right, but who has time to keep learning advanced calculus after college? (Well, I still study math & physics for fun, but over the course of decades, not years.) Not being able to see the world through these lenses I think means missing key engineering perspectives and relationships. Anyway, very well written article.

lordnacho 5y ago
I feel the same. How can it be that I went to a world famous institution providing 2-to-1 student-teacher ratios, but I still think the best explanations are these modern internet explanations? I guess the best explanations just bubble up in the modern environment. > you're not supposed to learn everything in college, you're supposed to learn how to learn But to learn how to learn, you gotta learn some things to a somewhat decent degree. I think at some point you need to have these linalg/divgradcurl things down, if only briefly. You might forget any particular topic, but if you've indexed it you should be able to pick it up again, particularly in the modern learning environment. Just imagine coding without access to StackOverflow.
- pfortuny 5y ago
  Honest question. 2-to-1 student-teacher ratio?? Where is that? Totally stunned.
  - ectopod 5y ago
    Oxford tutorials do this. They typically have one, two or three students.

panic_on_oops 5y ago

Wonderful post, thank you OP!

tobmlt 5y ago

See also from Keenan Crane and company (discrete differential geometry): http://ddg.cs.columbia.edu/SGP2014/LaplaceBeltrami.pdf http://ddg.cs.columbia.edu/SGP2014/LaplaceBeltrami.pdf “The Swiss Army knife of geometric operators.” I always thought that was cool since I usually think of diffusion in the context of fluid flow.

niffydroid 5y ago

You lost me at equation

Chris2048 5y ago

Reminds me of 3Blue1Brown "Divergence and curl": https://www.youtube.com/watch?v=rB83DpBJQsE https://www.youtube.com/watch?v=rB83DpBJQsE

s-macke 5y ago

You can also use it to solve labyrinths. Just put a high pressure at the beginning and a low pressure at the end. Solve the Poisson equation. The path through the labyrinth is always the steepest slope. In [1] you can see a small implementation of the idea. [1] https://simulationcorner.net/maze/ https://simulationcorner.net/maze/

setr 5y ago
More generally, I believe this is exactly the same as flow-field pathfinding, which conveniently has the property that you end up with a solution to the whole map -- so you've solved all pathfinding for all entities anywhere in the map, in one shot. Great for hosting a ridiculous amount of entities on a static (or slow-changing) map. Brogue's creator also "invented" djikstra maps, which I believe is also exactly the same[0][1], to handle AI strategic pathfinding (e.g. avoid hazards while reaching treasure). In Game AI Pro (I forget which article/book), someone suggests taking AI preferences (e.g. hunger, health) and computing a flow-field per preference (e.g. a food-map, a danger-map, etc), and multiplying the values against the preference to act as a weight (so food-map * %hunger, danger-map * %health), and summing the maps together to produce the final map used for pathfinding. Notably, the food-map can be shared by all entities consuming the same kind of food -- only the weight has to be re-calculated, and the final sum. [0] http://www.roguebasin.com/index.php?title=The_Incredible_Power_of_Dijkstra_Maps http://www.roguebasin.com/index.php?title=The_Incredible_Pow... [1] http://www.roguebasin.com/index.php/Dijkstra_Maps_Visualized http://www.roguebasin.com/index.php/Dijkstra_Maps_Visualized

alisonkisk 5y ago

The article never explains why ∇2 means "average of neighbors". It's the divergence of the gradient, which is (one kind of) n-dimensional 2nd derivative. In a one-dimensional function, the second derivative is 0 when there is no curvature, aka a straight line (of any slope), and any point on a line is equal to the average of its neighborhood. A plane also has this property, but in 2+ dimensions you can also make other shapes (like saddles), that are made up of lines (like a plane) but the lines are twisted relative to each other in interesting ways (like "string art"). You can also visualize (aka impose a coordinate system for) these surfaces as having positive curvature (concave up) in one direction, and exactly opposite negative curvature (convex up, or concave down) in the orthogonal direction, summing to 0.

vlmutolo 5y ago

> It is customary when simulating heat flow to use a wacky color palette where red is hot and blue is cold, with all kinds of intermediate colors in between In an otherwise excellent article, this is the only issue I could find. We really need to stop using/recommending/normalizing rainbow color maps (i.e. jet). They actively confuse readers by creating visual artifacts that aren't actually in the data. This article has some great explanations and visuals. https://jakevdp.github.io/blog/2014/10/16/how-bad-is-your-colormap/ https://jakevdp.github.io/blog/2014/10/16/how-bad-is-your-co... The original post uses a rainbow color map to represent temperature-related things because having a diverging color map is often a useful intuition for temperature heat maps. But in that case, we should prefer one of the following diverging color maps listed on the matplotlib site (this list definitely isn't exhaustive, but it is helpful). https://matplotlib.org/stable/_images/sphx_glr_colormaps_004.png https://matplotlib.org/stable/_images/sphx_glr_colormaps_004... More on matplotlib's well-chosen color maps: https://matplotlib.org/stable/tutorials/colors/colormaps.html https://matplotlib.org/stable/tutorials/colors/colormaps.htm...

just_temp 5y ago
Can not agree with this more, if people want to plot something that is linear please use a perceptually linear colormap! Just a one second glance at the Mona Lisa in rainbow/Jet is enough to make you gouge your eyes out. https://peterjamesthomas.com/2017/09/15/hurricanes-and-data-visualisation-part-ib-the-mona-lisa/ https://peterjamesthomas.com/2017/09/15/hurricanes-and-data-... For a more technical description the information behind the newer matplotlib defaults, particularly the scipy talk, is great. https://bids.github.io/colormap/ https://bids.github.io/colormap/ And for those that do not like the matplotlib options, colorcet provides a wider range of alternatives that are not trash (unlike Jet) https://colorcet.holoviz.org/index.html https://colorcet.holoviz.org/index.html
- vlmutolo 5y ago
  Thanks for the links! I especially like the Mona Lisa example. I'll probably steal that the next time this topic comes up.

cornel_io 5y ago

> This oval (called a separatrix) has the special property that no wind flows through it at all. It acts just like a solid surface. We have already assembled a reasonably accurate simulation of how air flows around an ellipse in some confined space like a wind tunnel! My understanding of airflow simulation (undergrad-level at best) is that the correct boundary condition is almost without exception the no-slip one: air should be stationary at the surface of each object, not just have zero flow across it. Am I correct that the calculation mentioned above really only applies to the "dry liquid" scenario where there is no drag and zero viscosity?

mferraro89 5y ago
You are correct! The simulation as written does not handle the no slip condition. The simulated "air" is inviscid and irrotational. You would need to tackle a few more things in order to have a really accurate simulation.
- SubiculumCode 5y ago
  Thanks for the awesome write-up.writer's note. The first paragraph of your conclusion would have served better as your introductory paragraph.

sgarrity 5y ago

I thought poisson's equation was expressed as: <><

toolslive 5y ago

I consider Pagerank to be a discrete variation of the poission equation.

Synaesthesia 5y ago

We referred to this method of numerically solving Poisson's equation by successively averaging values, as the method of relaxation.

Lichtso 5y ago

Last year a Monte Carlo approach to estimate the solution quickly was discovered. It works somewhat similar to diffusion curves. https://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryProcessing/index.html https://www.cs.cmu.edu/~kmcrane/Projects/MonteCarloGeometryP...

WalterBright 5y ago

A very understandable explanation of it, and its uses.

Robotbeat 5y ago

I have used Poisson’s Equation recently. Actually, something similar. Been coding up some thermal conductivity calculations form scratch. What’s nice about the Laplace equation is that there are some exact analytical solutions in some situations which is really helpful for validating numerical codes and sanity checking. It’s also simple enough to implement in an Excel spreadsheet, Color coding the cells to indicate temperature. Good sort of “no code” example of the concept.

daleroberts 5y ago

Here is my implementation in python of the Poisson equation on an arbitrary 2D domain using the finite element method. I used this for teaching a course in partial differential equations: https://github.com/daleroberts/poisson https://github.com/daleroberts/poisson

MauranKilom 5y ago

> From here we could use Bernoulli's equation to find the pressure distribution, which we could integrate over the surface to find drag and lift and so on. With a few tweaks we could simulate rotational flow, vortex panels, real wing profiles, and so on. > With just a few simple building block we're already edging up on real computational fluid dynamics. All this just by adding up some matrices! Correct me if I'm wrong, but the only "real" CFD you could solve with this are incompressible potential flows [0]. Solving Navier Stokes is clearly not just "a few tweaks away" from the Laplace equation, but I would be curious which tweaks would take you to e.g. rotational flows. [0]: https://en.wikipedia.org/wiki/Potential_flow https://en.wikipedia.org/wiki/Potential_flow

powderpig 5y ago

I too hate the coloured plots but I can tell you as a Stress & Structures engineer who works with ANSYS regularly, it does help when assessing material limits and strain energies. This is especially true when you're working with models where temperature is time dependent. A great article, one to bookmark for sure.

woopwoop 5y ago

Why is the laplacian so ubiquitous? Well, locally any reasonable PDE is well approximated by a linear one. And linear PDEs are of the form Lu = f for a linear differentiable operator L. But why does the particular case of L = laplacian show up so often in physics? Galilean invariance says that the laws of physics should be invariant under translation and rotation. Checking this against Lu = f, you can verify that one requires that L commute with translations and rotations. That is L(u_h) = (Lu)_h, where u_h(x) = u(x+h), and similarly for rotations. What you can show (pretty easily on the Fourier side) is that every linear differential operator with these properties is a polynomial in the Laplacian.

amirkdv 5y ago
> you can show [...] every linear differential operator [that commutes with translations and rotations] is a polynomial in the Laplacian. I may be tripping over something here but this doesn't sound right. If you mean polynomial in the real number sense, i.e. Lu = a0 + a1 Delta(u) + a2 Delta(u)^2 + ... (where Delta = Laplacian, and a0, a1, ... real numbers), then is this true? The famous wave operator doesn't have this form. And if you mean "polynomial" as in a series over function space, i.e. Lu ~ a0 + a1 Du + a2 D^2u + ... (infinite terms, not equality but convergence) where D is the usual differential operator and a0 is a number, a1 is a 1-d vector, a2 is a 2d matrix, so on, then this is standard calculus of variations on any reasonable function space. Taylor series for function spaces if you want. I don't think that's limited to nice Galilean operators.
- woopwoop 5y ago
  Thanks for pointing this out, I should have been more precise. I was referring to time invariant equations. You are right that this doesn't apply to the wave operator (which also does not commute with rotations in R^4). But in space only, yes polynomial means polynomial with constant coefficients. The proof basically goes that commuting with translations implies immediately that L has constant coefficients. Then on the Fourier side applying L translates to multiplying by a polynomial, and commuting with rotations translates to the claim that that polynomial is rotation invariant. And every rotation invariant polynomial in several variables is p(|x|^2) for some polynomial p in on variable. Then on the spacial side p(|x|^2) translates to L = p(laplacian).

Attained 5y ago

What does it mean? I fully expected to explain the semantics and yet I still don't know what that triangle squared means. Oh well.

steve76 5y ago

This article details matrices as fixed known sets of data. Much more interesting is the vector calculus. Vectors as the result of the sum of their components are good. But think of the reverse. Your model outputs the current vector. Your model predicts the next vector in sequence. Just from current and next, your model provides any amount of underlying components, as high resolution as you want. Nature fails before your model does, through the lack of continuity in steps and impulses and the mean value theorem no longer holds. Vector analysis was pioneered by Josiah Gibbs, someone who doesn't nearly get enough credit. Up there with Boltzmann, Maxwell, Planck, and Einstein.

carodgers 5y ago

Beautifully done. Everything was clear and intelligible. I have a CS background, but no physics/aerodynamics background. What are the minimum additional steps beyond this tutorial which could produce an aerodynamically correct 2D wing simulator? (With turbulence, not a steady-state solution.) I've come back to tackle this topic intermittently but have never cracked it. To anyone with expertise here who could share an overview and links, I'd be grateful.

btrettel 5y ago
There are many different ways to do what you'd like. The easiest starting point would probably be this tutorial: https://github.com/barbagroup/CFDPython https://github.com/barbagroup/CFDPython But that won't handle turbulence. The real "turbulence problem" is that computing actual turbulent flows requires enormous computational resources. So instead of solving the Navier-Stokes equations, related equations with lower computational cost are solved. Because of how these equations are developed, they require modeling of "unclosed" terms, and this is a likely source of inaccuracy. If you want something relatively simple, you could take the RANS approach and use the Spalart-Allmaras model: https://www.cfd-online.com/Wiki/Introduction_to_turbulence/Reynolds_averaged_equations https://www.cfd-online.com/Wiki/Introduction_to_turbulence/R... https://www.cfd-online.com/Wiki/Spalart-Allmaras_model https://www.cfd-online.com/Wiki/Spalart-Allmaras_model How to implement the changes to the final part of Lorena Barba's tutorial should be fairly obvious by the time you get there.
- carodgers 5y ago
  The CFD Online site is a great resource that I had missed before. Thank you for the links.

ybarhighbar 5y ago

Always looking to interrelate my knowledge that sometimes appears to more disparate that I care to consider. The simplicity of considering a matrix function in this way is very elegant. This type of methodology can help demystify some of the earlier steps in modeling by adding a simple intuitive picture that encourages a connection to calculation from the get go.

aabbcc1241 5y ago

The article explain clearly. Made a toy canvas to play with poisson https://github.com/beenotung/poisson-canvas https://github.com/beenotung/poisson-canvas

amai 5y ago

"you should know that Laplace's Equation has a close cousin called the Biharmonic Equation... Which is widely used in modeling quantum mechanics. " Can someone elaborate on that? The Wikipedia article (https://en.wikipedia.org/wiki/Biharmonic_equation https://en.wikipedia.org/wiki/Biharmonic_equation) does not mention anything about quantum mechanics use cases for the biharmonic equation.

vayhay 5y ago

Hiện tại hệ thống VayHay đang tuyển nhân viên Telesale khu vực cần thơ và các tỉnh thành khác toàn quốc làm việc 8h/ngày và có hỗ trợ vay tiêu dùng toàn quốc không thế chấp kì hạn và các gói vay linh hoạt, truy cập VayHay.vn để biết thêm thông tin