top | item 27748967

Poisson's Equation

766 points| mferraro89 | 4 years ago |mattferraro.dev

159 comments

order

jbay808|4 years 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-phot...

codethief|4 years 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://minimalsurfaces.blog/ (lots of illustrations)

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

mixedmath|4 years 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.

woopwoop|4 years 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|4 years 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.

prionassembly|4 years 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/

throwamon|4 years 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?

chris_st|4 years 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...

alisonkisk|4 years 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.

pphysch|4 years 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.

zitterbewegung|4 years ago

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

hasmanean|4 years ago

Apropos of nothing, some of my first programming lessons were FORTRAN programs from my dads old engineering textbooks. The book’s graphics showed them being handwritten on graph paper and even showed the punch cards they would ultimately be put on.

One of the most impressive programs solved Laplaces equation for heat flow in a pipe. They used some mysterious plotting library that put ASCII art to draw contour lines and output directly to a line printer. I thought it was the coolest thing ever, but it was hard for a high school kid to understand.

Over the years I’ve thought of writing that program in a modern language, but it was never worth the time and effort to go over that FORTRAN program in detail. Numerical methods are too often explained from the point of view of mathematicians and engineers, and not computer programmers.

The blog post here is 1000x more readable and much higher quality and does a lot to demystify the subject. I especially like his progression from simple brute force methods to more efficient solutions…which is the natural way to learn. His comment about “Laplace’s equation just means every point is the average of its neighbours” is perfect.

Believe or it not I was in need of just such a tutorial … for something I’ve been thinking about at work. His article really hit the spot and I’m grateful for it.

cpp_frog|4 years ago

It is easier to study from a theoretical point of view (easier that the heat or the wave equation: [1], Ch. 2), and it's easier to implement a solving method. When you are learning the finite element method, this is one of the first examples that people use to test if they got the right implementation.

Now, I wonder if the author regards it important in his particular area (aerospace engineering), I'm new to the field so I don't see how. Right now I'm reading a book [0] on models to solve problems concerning aerospace applications and they mostly use a simplified form of the Navier-Stokes equations together with some elasticity assumptions.

[0] Fluid Structure Interacion, Morand-Ohayon.

[1] Partial Differential Equations, Evans.

nyghtly|4 years ago

The opening of this article had me recalling an old TED talk: "People don't buy what you do, they buy why you do it."

Of course, the author here isn't really selling anything, at least not to the lay man. He's writing niche articles for a niche audience. That is to say, not me.

https://www.ted.com/talks/simon_sinek_how_great_leaders_insp...

OnlyOneCannolo|4 years ago

It's a simple and general method that works in many domains, which is an usual combination.

The other ways of solving the example of arbitrary heat sources and sinks on a plate range from hacky combinations of simpler methods, to tedious math, to complicated general methods. If you switch from heat to pressure distribution, you'd have the same types of options, but the specific methods would be different.

wildmanx|4 years 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.

s-macke|4 years 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/

setr|4 years 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_Pow...

[1] http://www.roguebasin.com/index.php/Dijkstra_Maps_Visualized

vlmutolo|4 years 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-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...

More on matplotlib's well-chosen color maps:

https://matplotlib.org/stable/tutorials/colors/colormaps.htm...

just_temp|4 years 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-...

For a more technical description the information behind the newer matplotlib defaults, particularly the scipy talk, is great. 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

daleroberts|4 years 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

Robotbeat|4 years 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.

enriquto|4 years 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.

yngvizzle|4 years 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.

phkahler|4 years 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.

Porygon|4 years ago

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

pgustafs|4 years 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|4 years 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.

MauranKilom|4 years 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

SavantIdiot|4 years 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|4 years 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.

amoorthy|4 years ago

Agreed! It's one reason I'm bullish about our children learning and applying math better than us.

Btw, 3 Blue 1 Brown's videos on linear algebra [1] are similarly awesome and of course his video on Fourier is magnificent [2]. Another awesome math explanation is on game theory by Nicky case [3].

[1] https://www.youtube.com/playlist?list=PLZHQObOWTQDPD3MizzM2x.... (you'll finally learn what a determinant is!) [2] https://www.youtube.com/watch?v=spUNpyF58BY [3] https://ncase.me/trust/

senthil_rajasek|4 years ago

30 yrs ago the internet was still in its infancy and good engineering teachers were available to only a chosen few and closed out to the rest of the world.

I stayed curious and learning. Better and better teaching methods became accessible thanks to the maturing internet. I started to understand ideas in engineering much better.

Alas, there is no way I can express this progress in a resume but learning and understanding satisfies my curiosity.

p_j_w|4 years ago

>I could have used this 32(!) years ago when I was struggling in college.

I was struggling with this material more like 15 years ago, but same. I wish someone had explained the LaPlacian like this when I was in Multivariable calculus:

>Find me a function f where every value everywhere is the average of the values around it.

It's so simple and easy to grasp, yet provides so much insight into what's going on when you're doing the actual calculation behind the operation, but is so easy to lose sight of when you're overwhelmed with figuring out the 'mechanics' of it and everything else that was covered in the day's lesson plan.

MontyCarloHall|4 years ago

> rather than add a curve to a test that makes a 23 out of 100 an "A" grade.

For people wondering why fresh college grads they interview somehow have 4.0 averages yet can’t code FizzBuzz during an interview, this is the answer.

It’s also why good hiring managers do not even bother to look at GPAs listed on CVs. They are so inflated as to be totally meaningless.

bradstewart|4 years ago

For me at least, at a large public university with a highly-rated engineering department, it mostly came down to math courses being taught by TAs/new faculty who were not primarily interested in teaching.

When I had actual engineering dept courses on these topics (frequency domain comes to mind), I grok'ed it pretty quickly. But the math lectures were basically spent doing practice problems from the book, without any sort of insight or practical application.

I just wanted to build shit in college.

jozvolskyef|4 years ago

The generation growing up today is the first that can broadly learn mathematics from mathematicians rather than from teachers, and we should encourage them to do so.

When I was 12-16, I was interested in physics enough that I would study it in my free time just for fun, solving problems from the Physics Olympiad. I could solve lots of the problems with just intuition, until I stumbled upon the sliding chain problem. I spent about a week pondering over it with no results. I approached my math teacher. He admitted defeat and referred me to the physics teacher. If I remember correctly, the physics teacher avoided admitting he didn't know how to solve it and didn't give me any pointers. That was the last Physics Olympiad problem I tried to solve, after dreaming of attending the competition for years (I should mention that at that time my primary interest already were computers; physics was just a hobby). I wonder whether this could still be a problem today when anything can be found on YouTube. There's a lot of noise, too. How long does it take the average curious kid to find the signal that is 3b1b, MIT OCW, etc?

matheusmoreira|4 years ago

> 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

School is not about learning. It's about money. Seating as many students as possible for as long as possible to suck in as much student loan money as possible. For this you need to be a prestigious school. Of course they want to boost their grades.

These are the people who have the gall to claim a moral high ground when they find and punish "cheaters".

H8crilA|4 years ago

What's a "test curve"? How does it make 23/100 be an A grade, is it some kind of CDF-based transformation that makes certain % of students pass the test?

alisonkisk|4 years ago

Your 30 years of experience is why the article makes so much sense. If you saw this in college you wouldn't have also learned everything else about div and grad.

> these concepts are learned as succinctly

? Succinctness is *why" you didn't learn multiple interpretations of everything.

cornel_io|4 years 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|4 years 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.

amai|4 years 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) does not mention anything about quantum mechanics use cases for the biharmonic equation.

pm90|4 years ago

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

davidkuhta|4 years ago

> probably

can lead you to some fun places.

sgarrity|4 years ago

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

carodgers|4 years 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|4 years 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

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/R...

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.

ybarhighbar|4 years 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.

Attained|4 years 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.

powderpig|4 years 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.

toolslive|4 years ago

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

Synaesthesia|4 years ago

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

seemslegit|4 years ago

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

WalterBright|4 years ago

A very understandable explanation of it, and its uses.

vayhay|4 years ago

[deleted]

Jeff_Brown|4 years 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|4 years 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|4 years ago

A matrix can be seen as the discrete representation of a function…

jungturk|4 years 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.

gmmeyer|4 years ago

Poisson's Equation takes another function as an input, the matrix is the representation of the output of said function