Good luck with that. “Old Man River” has been avulsing and winding wherever he pleases for a few tens of millions of years, and trying to dictate its course is like pissing into the wind. John McPhee had a few good words about that as I recall.
Stranger
I am sure that someone, somewhere is trying to invert an ill-conditioned matrix.
That, and any sort of work done near a coordinate singularity (or a real singularity, of course, but those are much rarer). I recall one fun project I did in an orbital mechanics class: We were given a list of orbital elements with error bars on them, plus a list of definite measurements (with assumed zero error, but not enough of them to fully determine an orbit on their own), and had to use those to produce a better estimate of the true orbital elements. The catch was that the eccentricity we were given had an error estimate larger than the eccentricity itself (i.e., the error bars extended into the negatives), even though eccentricity is by definition always positive. Most students’ codes blew up pretty badly on that.
Later, in grad school, I had a numerical methods class, one of those fun ones where the professor was figuring it out alongside the students. The textbook talked about a certain class of differential equations, “stiff” equations, which would resist solving by most conventional techniques, and included special techniques that would work for them. But none of the examples in the textbook were actually stiff, leading all of us (including the professor) a bit puzzled about just when those techniques were actually called for. The next unit, the prof gave us an exercise based on an actual physical situation, IIRC something involving convection cells in the atmosphere, and none of us were able to solve it. Turns out, that was actually a stiff equation.
The worst operation is subtraction. The error is additive but the caculation result is deminished. I recall one example - forget how he got there - where the prof shows the result in a quadratic equation, and the error in the first term is such that as he said “the error is greater that the coefficient so you don’t know if this describes a parabola opening upwards or down.” The same applied to calulating the intersection of nearly-parallel lines.
A common issue for this sort of error growth is a repetitive loop, where the value is being refined without significant growth yet the error is additive, adding more error as the loop iterates. This is a problem with computers, which can happily do dozens or hundreds of iterations of a loop, something which humans rarely do by hand.
This is part of what @Dr.Strangelove was referring to when he said that it was a complicated subject: Some iterative loops will grow error, and some will shrink it. In fact, for a loop of one type, there’s often a way to rearrange the calculations to turn it into the other type. For instance, for the golden ration \phi \approx 1.61803, it’s true that \frac{1}{\phi}+1 = \phi. If you start with any random number, and alternate taking the reciprocal and adding 1, you’ll fairly quickly converge to \phi. But by the same token, it’s also true that \frac{1}{\phi-1}=\phi. If you start with exactly \phi, and repeatedly subtract 1 and then take the reciprocal, you’ll stay on exactly \phi, but if you’re even a little bit off (say, from finite floating point precision), the error will grow until you go completely crazy.
It’s just another example of humans ignoring nature, kind of like building on a shoreline frequently hit by hurricanes or floods, or building on a cliff. They have managed to keep the river more or less in its banks - but it’s expensive.
I lived in Louisiana for a while and they are well aware of the problem.
That’s not too common in my line of work, where most stuff is done in boring low-dimensional Cartesian coordinates with matrix math (where maybe at the end we do a one-time projection from 4D down to 3D). The one exception would be Euler angles, and the inevitable answer is: for the love of God, don’t use those. Use quaternions instead. You don’t have to understand them; just find a library. So weird coordinate singularities don’t play a huge role.
Excellent example. That makes me wonder if there’s some kind of symmetry property, or mathematical dualism, where you can always convert an iteration with shrinking error into growing error.
(Also, I see that I’m not the only one that regularly types “ration” for the word “ratio”. I know the dang word, but my fingers love that extra N).
Well, the glib answer would be “Invert your step, so whatever you’re using to get from x_n to x_{n+1} instead takes you from x_{n+1} to x_n”. But of course, not all mathematical operations are easily inverted.
I should also addendate my previous post: Starting near \phi and repeatedly subtracting and inverting will, in fact, eventually go crazy… but just after it goes crazy, it’ll start converging again to the other attractor, at 1-\phi.
(and I didn’t even realize that I was prone to that typo. Or maybe I’m not, and just did it that one time. I’ll have to pay attention now.)
I tried it for the common iteration x' = x/2 + 1/x for \sqrt{2}, but I just got x^2=2, which wasn’t very interesting. Maybe it only works when the variable shows up once on both sides of the equation.
I remembered one moderately common formula that has numeric stability problems: Heron’s formula. That is, to compute the area of a triangle from its sides:
s=\frac{1}{2}(a+b+c)
A=\sqrt{s(s-a)(s-b)(s-c)}
All those subtractions are killer for long, thin triangles. If a=b=1000000000000 and c=0.3, say, you can lose most or all of your significant bits.
You didn’t invert it properly. The inversion of that would be x' = x\pm\sqrt{x^2-2} (plus or minus depending on whether you’re starting off above or below sqrt(2)). Which does indeed grow out of control, very quickly and very badly.
EDIT: I should also mention, of course, that if you’re capable of doing that iteration, then you obviously already have some way of doing square roots, so the practical thing to do would be to just calculate sqrt(2) using that method directly.
I’ve forgotten the details, but when I first started using C, which didn’t support complex numbers, I implemented them manually. If you naively implement complex division, you can easily get catastrophic cancellation under some circumstances. It is easy to fix, but it took me a while to figure out what I was doing wrong.
Hmm… (13+21i) / (21+34i) should be very close to 1/phi+0i. Let’s see what happens.
\frac{13+21i}{21+34i}
\frac{13+21i}{21+34i}*\frac{21-34i}{21-34i}
\frac{273-442i+441i+714}{441+1156}
\frac{987-i}{1597}
0.6180338134…-i/1597
Yeah, the imaginary parts almost cancel out in step 3, but that’s not so big a deal, since the result of that is just a small number, which is what it’s supposed to be.
With some more thought: The place where near-total cancellation would be a real problem would be in a denominator. But since the denominator is the product of a number and its complex conjugate, it’ll always be a sum of two squares, not a difference. So that shouldn’t cause major problems. The only way it’d ever be small would be if your original denominator was, and small denominators causing problems isn’t anything special to complex numbers.
I’m curious what special cases you found that caused the problems.
It’s possible to compute a complex square root of (a+bi) as:
u = \pm \sqrt{\frac{a+\sqrt{a^2+b^2}}{2}}
root = (u + \frac{b}{2u}i)
If a is large and negative and b is small, you can get catastrophic cancellation. The answer is… don’t do that (use De Moivre’s Theorem, adjust the formula based on the size of the numbers, etc.)
But people should understand them. They’re incredibly useful and turn really complex vector math in spherical trigonometry into simple arithmetic and basic linear algebra.
Stranger
I’ve been fascinated with them since I discovered them in high school, and loved them even more when I learned how elegantly they could represent 3D rotations. So I certainly agree in that sense. Sadly, not everyone seems to have the same fascination with them, or with any math for that matter. So to them I say to pick up a library. They’ll work even if you don’t know how they function on the inside.
Sometimes Euler angles are a part of an actual physical system, though. Take a telescope mount, for instance. Most telescope mounts have two pivots, at right angles to each other (either azimuth and altitude, or right ascension and declination). If you were to add a third degree of freedom, to roll the tube of the telescope around its axis (say, to determine the orientation of an image on the sensor), you’d now have a 3D object, the telescope, that could be put into any orientation in 3D space. And those three degrees of freedom you had would, in fact, be Euler angles. Which does lead to difficulties, when you’re working near one of the poles (either the celestial poles, for an equatorial mount, or the zenith, for an alt-az).
I’m sure you could build a telescope mount based on quaternions, instead, but it’d be mechanically much more complicated.
Well, that’s true. But it’s also true that the coordinate singularity is also a “here’s where your telescope breaks” point.
It would be difficult using analog dials for measurement and adjustment but really easy using a software based orientation system on an unconstrained ball joint. Most modern manipulative robots and essentially all inertial orientation systems like those in hour smartphone use quaternions to calculate orientations and articulation mechanics because trying to handle something with completely arbitrary dynamics or a system with more than six unconstrained degrees of freedom would be confusing and error prone. And we hate, hate, hate mathematical ‘gimbal lock’ when there is no actual physical constraint which produces it.
Stranger