rambling

I tried to avoid continuous case, but yes. The real heat transfer is great. ∂T/∂t = c ∂^2T/∂x^2 (change of the temperature in time is equal to the second derivative in space*)). From theoretical standpoint everything works, you hit it with Fourier transform (that equation is called Fourier equation, and this is not exactly a coincidence) some of the strange symbols disappear and we see that the evolution of heat distribution is just smoothing it out. As we sould expect from heat. euler scheme
But if we want to calculate if for real, we divide time into segments to look only at discrete points in time, and divide space. For example we look at discrete points (the state is represented by values in those points), and the strange letters become differences,  ∂T/∂x -> 1/dx(T[n]-T[n-1])  ∂^2T/∂x^2 -> 1/dx(T[n+1]-2T[n]+T[n-1])  (at least in 1D, 2D is a bit longer). FEM - finite difference method. There are better (sort of) methods**) but in a regular square grid and low order methods the equations turns out the same.
Now, both those discretization will work together to beat us in the shiny metal...
I have already written the space part with unnecessary details, I will leave it as a bonus for masochists. The important part it, taking all the temperatures we calculate the change of temperture. dT/dt = f (T). That function, if we forget about the "difference >1K for transfer") is linear. So it is a matrix dT/dt = A * T.
Now, the main villan, discretization in time.  They most likley use some _explict_ scheme. It just means we can calculate the new values of T using values of T from previous iterations. At first it depended on the order of pipes, but not it does not, so it is explic euler, maybe applied in chessboard  pattern. Lets say it is just euler. If we have chenge (rato of the chenge) dT/dt = A T and the initial state is just T, whet is the state after dt time? A good guess is
T_new = T_old  + A T_old. This is called Euler explict scheme

Now the whole eingenvalues/eingenvectors stuff come in handy if vk is an eingenvector with e.value ak, then A vk = ak vk. So, is we apply the evolution to heat distributed like in eingenvector, the change in time is proportional to the same vector. It gwows or shrinks, but does not chenge shape!
We can analyse the evolution not by looking at each temperature, but by looking at the teble of temperatures as a combination of all those eingenvector. The vector T = sum_k bk vk  where bk are just some real coeficient, telling us how much of a given vector is in the miv. And since the evolution of each one is saparated, we just get formula on how the numbers bk chenge.
Apply our expression for T into the equation, separate different vectors and we get
bk_new = bk_old + ak bk_old = bk_old (1+ak)
In the continous case we also can use eingenvectors. All (but one that represent the flat steady state) have negative eingenvalues. This mean the coeficient expotentially vanish. We are left with the steady state, all pipes have the same temperature (we did not put sources yet, but with sources it similarly gets to a steady state). Will the discrete version do the same?

Lets skip forward n ticks. Then bk_n = bk_0 (1+ak)^n.
 We have said the eingenvalues are negative. If they are negative and small everything is great. (1+ak) is smaller than 1, but greater than 0. After each tick the coeficient bk is slightly smaller. This is what we want.
But what if ak is too negative? In the continous case it just mean that mode it dissapears very quickly.

But in discrete case, if  ak is smaller than -1 (but greater than -2)
 (1+ak) is negative. The amplitude is decreasing, but the temperatures oscilate. It looks wrong.
If ak < -2, then (ak+1) <-1. This not only oscilates, but the abs(bk) goes biger and bigger. It is catastrophic for our simulation.
OK, so what the ingenvalue ak depend on. First, heat transfer coeficint. Smaller coeficinet, smaller chenges. Second, the time step. It is like dt. Smaller step size means smaller changes in that time. And, for last, the eingenvalue of laplacian. For a line of pipes it gets close to -4/h^2. his the size of the space grid, we can put it h=1 (the heat transfer  coefficient has to be expressed in the same units for this to work). We and up with

ak = -4*coeff*dt

An, and the coeficinet is expressed in ticks also, so we can get rid of dt (dt=1)
ak = -4 * coeff
In the previous post I linked a discusion where devs said the coef was chenged from 1 to 1/5. -4/5 = -0.8
(1+(-0.8)) = 0.2 No oscilations, dissapering modes, great success;-)
The only problem was, for a 2d grid that most negative eigenvalue is two times greater (from memory). We went back to oscillating, but still damped regime.

I have to stress it out. The above is a tool for analysis. To see what will happen. Not how it is computed in the code. There is no eingenvectors used, all calculations are directly on the table with heatpipes and temperatures. Most likley there is no sparse matrix expressed directly, just information who touch who.

No, this is not about floating point precision. They get bad rep constantly (pipes in satisfactory was off, floating point was a scapegoat too).
But yes, you may see it as too much heat is transfered per tick. In the continous case we have a nice exp(-at) curve. But now we are simulating it in discrete time. Choose a point on the curve, take a tangent line to it, and move along the line by dt time. If the dt is too long, we escape to the other side of the graph (below 0).

How to deal with it? In real similation, you may apply implic numerical schemes. They try to use the slope of the function from the "now" point. Implic euler is always stable. But we do not know where that point is, so we and up with a system of equation. In a bare form unacceptable in game due to computational cost. It is sparse system (if the player is sane), so may be still O(n) and we can get away with calculating it less often then every tick... it is a big maybe. I would rather go with tricks, the chagnges are slow (replace a long pipe with a one element with enough states so it looks ok) or use the fact, we are probabli limited by access to cold memory, so a but more computetion can be done, as far as we do not read/write too much.
Also, physical parameters can be tinkered with. The ak eingenvalue depend on the speed of transportation and heat capability. But we do not want to decrease the first, nor probably to increase the second one (just like with pipes, they could not increase throughout significally without making pipes hold more fluid. This limitation was seen probably more directly - pipes drieds up - but it is all conected).
On the other hand, maybe there is just no need to do any chenges heatpipes.



Lets take a long heatpipe. N segments. We have a table of temperatures, T[0].. T[N-1]. What is the change after a tick?
const * (T[n+1]-2T[n]+T[n-1])  ***)
If we write it as a matrix, we get somathing like this image (782×444) (d2vlcm61l7u1fs.cloudfront.net) - that one is on a circle, the last pipe touch the first, so there are 1 in the corners
(-2 in the diagonal, 1 above and below, 0 elsewhere). Lets call it A.
A is a very nice matrix. The eingenvalues are -4 sin(pi k /(2(n+1)))^2. The one with the biggest amplitude is near -4. The smallest around -pi/(n+1)^2


*) why second? Flux (power escaping in given direction) if proportional to temperature difference. And the change of energy is a difference between flux a bit to the left and a bit to the right.
**)Or we divide the space into shapes, most likely triangles, and we say temperature in those shapes can have its own shape (for example, it can be a plane). Then we see how that entire shape will play with the equation. skipping weak formulation and weak solutions...   This is the famous finite element method. Or Galerkin.
***) The code in factorio may do it in a bit different way. It doesn't matter here, the end result is similar. And almost for sure they do not use any matrices, just loop through the table of temperatures.