?

Log in

No account? Create an account
 
 
20 April 2008 @ 02:30 am
(Not) Only Solutions  
Okay, folks, this is another super-nerdy one; so, if you're not interested in math, this is your warning to scroll down.

Still with me? Good for you! Here's the thing: I love math, and while trying to learn how to simulate physics processes using JAVA, I come across a lot of problems that require it. You see, JAVA does simple calculations very well; but when a way is needed to transfer paper problems to JAVA calculations, thence to representations like two and three dimensional graphics, a lot of fudging has to be done. First off, JAVA doesn't come with matrix operations, although there are many packages in development for this, like JAMA. I'm a glutton for punishment, however, and insist on being able to do things on my own, so I have gone to the trouble to create classes for my own use.

Now, matrix math isn't remarkable by itself; it's a natural offshoot of (linear) algebra using two dimensional arrays rather than constants and vectors for variables. The differences between the two aren't trivial, however, and this is where a good text on abstract algebra comes in handy. I won't go into the details, but suffice it to say that, in matrix algebra, A * B isn't always equal to B * A, and 1/A may not exist, even if A is not the empty matrix (meaning that all elements aren't 0). This can get ugly when using matrices to solve linear equations, especially when doing so is required for otherwise straightforward physical problems like simulating nonlinear partial differential equations for a square "wave" tank. Apart from having to come up with array-based matrix math from scratch, I have to solve the general problems ahead on paper and use finite difference methods, because JAVA doesn't have built-in functions for symbolic calculus, either. In short, to write a graphical simulation of seemingly simple water waves (or the more complex yet still calculable patterns of light refracted through them), I have to reinvent the tools which others have done much more ably.

Programs like Mathematica (which I used extensively until it became too expensive to maintain) do all these things readily and easily, making graphical simulations practically child's play; but my ultimate aim isn't to create the simulations themselves, although the joy I get when finally getting one to work is sublime. I'm doing it to learn, like I do most things. Through this nonsense, I hope to get a more thorough understanding of programming, the relative advantages and disadvantages of particular strategies for simulating physical phenomena, and math itself. I have done as much for signal and image processing, building "GUI" desktop services, and trying my hand at developing a rudimentary symbolic math language like Mathematica; but those are epic stories in themselves, fraught with frustration and fleeting satisfaction. It has always been the underlying, esoteric math that has drawn me; and this is what I always come back to when other flights of fancy get too far afield. And here is where I came upon a pestering problem I've had with math ever since studying tensors and abstract algebra.

When manipulating vectors, that is, arrays of a single depth with several elements, like {A0, A1, A2} for a vector "A" of length 3, it's traditional and basically straightforward to use matrix math. One vector may be transformed by 'multiplying' it with a two dimensional matrix with the length squared number of elements, like {{B00, B01, B02}, {B10, B11, B12}, {B20, B21, B22}} for matrix "B" to create a new vector "C". This is done by adding parts of each that are multiplied in this way: B00 * A0 + B01 * A1 + B02 * A2 for the new vector element C0. For C1, we have B10 * A0 + B11 * A1 + B12 * A2; and for C2, B20 * A0 + B21 * A1 + B22 * A2. For brevity, this is normally written in the 'Einstein notation' or 'abstract index notation', B_n,m  A_m = C_n, where it's understood that the "B_n,m  A_m" stands for multiplication and the doubling of subscripts within a single term ("m" in this case) implies summing over the number of "C" elements. This is convenient in that it remains simple for vectors and matrices of any size (or tensors, spinors, or other exotic arrangements of elements). It's also convenient for me in that for a vector of a given length, each element of the matrix is used exactly once, so that for any transformation "A" to "C", there is only one corresponding matrix, if one exists at all. In a similar way, matrices are 'multiplied' with each other to create new matrices, but I won't go into that; except to say that applying two or more transformations to a given vector can be done by multiplying several matrices together first, then applying only the resulting one (and remembering that, for matrix multiplication, A * B rarely equals B * A).

Confused much? Laughing at my childish ramblings? Still interested? Let me show you how this kind of math might be applied to a simulation. Say, for instance, we are trying to show the motion of a string in tension, like that of a guitar or violin. For this problem, the 'one dimensional wave equation' is a ready approximation that can give exceedingly accurate results if solved correctly. The problem with this, in my case, is that it most readily gives continuous solutions, using sine and cosine functions for strings fixed at both ends; but computer simulations (especially those done with JAVA) become very slow and difficult to manage when the result is expected to be so accurate, and finding solutions by continuous methods can be extremely time consuming in itself (and you can forget generalizing them in the program for arbitrary initial or boundary conditions). What it usually done, and what I have recently been trying, is finite difference methods, in which the field (here, a simple length of string) is divided into regular lengths and the time is divided into finite fractions of a second. In this way, and for reasons that are too in depth to go into, all the calculations that are needed are simple additions and multiplications rather than higher order functions. In the case of a string that, for the foreseeable time, has no parts that overlap or become vertical, and is close to straight in reality, we can divide it into a number of segments of equal distance along the X axis and solve the given problem for each in turn. Where before we had a single, continuous curve to find, tempting though that sounds, now we need only approximate that curve with a number of joined line segments. If our computer can handle it, we can use a large number of such segments, even to having one for every pixel-length along the string shown. Remember, it's not necessarily the accuracy of the simulation that we're improving, but the amount and complexity of computations (because, in calculating higher functions, computers basically use the same finite methods, but with far more steps and to a much higher degree of accuracy than is needed for the graphics). All this means that, instead of finding a certain overall function for the curve, like Y = Cos(x) Cos(t) - 1/3 Cos(3 x) Cos(3 t) + 1/9 Cos(9 x) Cos(9 t)... and calculating each point along the length with that for every frame, we set up a 'vector' with our number of line segment end points {X0, X1, X2... } and calculate only those for the curve height values {Y0, Y1, Y2... } at given instants of time {T0, T1, T2... }. Using the fact that each X_n is of uniform distance and each T_m is also of equal length, we can use a matrix transformation to incrementally change the shape of the string from instant to instant. Y[t + 1] = A * Y[t]; for the vector of Y values and a matrix "A" to be found by solving the wave equation.

This is pretty rudimentary in programming and so I'm only doing a lot of work that I really don't have to-it's all been done for me. I'm not in it for the final product, though, as I said. I'm in it for learning how it's done. And I'm doing quite well with knocking out a few Rube Goldberg-esque programs along these lines. I have, in fact, written up ones for John Conway's Game of Life (with many variations), populations studies, galaxy simulations, fractals (of many types), and even billiards tables. How workable they are for JAVA applets doesn't really interest me except that I might put them up on a website for grits and shins someday. As long as they don't crash, and show something close enough to the 'real' thing for my taste, I'm happy. However, sometimes I come upon a curiosity that sends me on a wild tangent, and that's the reason for this post.

Calculating an approximation of a vibrating string is easy enough, but what about a wave tank? What about an oddly shaped one? And, more importantly, what about three dimensional volumes like sound wave simulations, or nonlinear systems like turbulent flows? Here, we would not be working with vectors of N points, but whole arrays of N^2 or N^3 points, not to mention higher dimensions for use with Relativistic systems and the like. For these types of things, the normal route is to again divide the 'field' into regular areas, a grid, and calculate the values for each point at a time. Here, though, there's a problem. For even a two dimensional grid on a square area, we are no longer using a vector of values, but a matrix! And for volumes of higher dimensions, what? A 'hyper' matrix? Matrix math only accommodates N * M rectangular arrays of elements, not N * M * K, etc. How are we supposed to calculate those? Well, usually, in those cases, the formulas for each point are calculated directly, without resorting to abstract algebra. For instance, Z[X_n, Y_n, T_m + 1] = -Z[X_n, Y_n, T_m - 1] + 2 (1 - (m/n)^2) Z[X_n, Y_n, T_m] +... , for each Z value on a square grid. I'm fine with doing it this way; it's straightforward enough and probably the best way, but the matrix method is more elegant. The fact that matrix multiplication can be used to transform over time the elements of a grid doesn't really help, though, because it doesn't transform the grid in the complete way that it would a vector. Let me give just a ghost of an explanation. in matrix-and-vector multiplication, every element of the transformed vector is created by a number of elements of the transform matrix equal to its length. That is, it's equivalent to every point along the length of the string having a relation to every other point, with that relation dictating its movement. That might sound obvious, but it's not. Another, badly posed, solution might have each element being used twice or more times in the calculation; the equivalent of a point in space whose movement is dictated by more than the number of points in that space! Or less. This is called an over-determined or under-determined system; and, though each can have its uses (in particular, wave function solutions actually come out as far under-determined, so their transformation matrices are sparse), it's not a true enough representation of the system for my uber-nerdy disposition. When matrices are multiplied with other matrices, as in a transformation, each element is only determined by an N-number of its neighbors, rather than all N * M members. In the notation; A_n,k B_k,m = C_n,m; or A00 * B00 + A01 * B10 + A02 * B20 = C00 for a 3x3 matrix. That's the equivalent of the motion of a point in a square only being determined by the points along a line.

That's not quite what it means, since it's only a simulation, but it's the loss of a large number of determining factors that is the problem. But my problem with it isn't why it isn't used for this. Equating for a square area with an infinite number of points (a continuous medium) couldn't be done this way, anyway. And there's no need for a more complex math when the straightforward approach works. Still, the first time I used matrices for calculating the elements of a string I automatically wondered how this might be done for higher dimensions, and considered matrices as well as tensors (though tensors are actually not for use with simple constant elements but are an entirely different thing). Even though I found a way to transform a matrix using two others, one as a 'left-hand operator' and the other 'right-hand', that would give a full N^2 transformation (the original matrix is transposed, or turned sideways, for the second multiplication), that's again fudging it just to make it fit. What I needed was hypermatrices, or N^(something other than two) arrays to use for transformations, but my texts are silent in that matter except for tensor algebra, which isn't the same thing. I needed a way to find the matrix-y things like inverses, Unit matrices, and determinants for these impossible matrices of three and higher dimensions. Well, I couldn't find any, so I invented them, as usual. I'm sure that it's been done before and thrown out as stupid, but I did it just the same.

It turns out that the simple operations like addition and scalar (constant) multiplication act with matrices of any dimension the same way, and a matrix of 2 * N dimensions can be multiplied with one of N dimensions to yield another one of N dimensions. In fact, this is the only way in which the transformation is 'complete'. An element of a square grid is determined by all the elements of that grid once when a matrix of four dimensions is used and an element of a cubic grid is equally perfectly determined when a matrix of six is used, etc. It was in trying to figure out how to do determinants and inverses that I came upon the elegant part that proved that it all would work and exactly with the same properties as normal matrix operations. If we 'condense' a 2 * N dimensional 'hyper' matrix into a regular two dimensional one by flattening out each N dimensional group of elements and do the same for its companion N dimensional 'hyper' vector, take whatever determinants, traces, and inverses etc. that are needed, then 'expand' them again, they transform perfectly! In other words, it's really all two dimensional matrix math, or is all based upon that. It only works for hypermatrices with twice the dimensions of their hypervectors, however, and determinants and inverses require 'hypercubic' hypermatrices just like only square matrices can have regular inverses. Yet, an irregularly shaped area can be approximated using a square matrix with special border elements and all zero elements outside the border, so any shaped space of any dimension can still be simulated using these hypermatrices.

The real problem, and this is the one that dooms my idea for any practical use, is that such complex beasts are so complex. A simulation of a string can be made to look very realistic using only a few, or a few dozen, elements, so that the matrix used to transform it also only uses a few, or a few hundred, elements. A two dimensional grid with a hundred elements on a side (tiny, even for use on a small computer screen), on the other hand, would require a hypermatrix with 100^4, or one hundred million elements! One needed to transform a three dimensional cube of the same size would need a trillion! My computer is fast, but it can't do tens of trillions of operations per second... yet. And using this method for a truly accurate simulation, as is used for weather or multiple body problems, would be absolutely beyond ridiculous. So, as far as my JAVA programming is concerned, hypermatrices are only a curiosity to be played with; but at least I know they work, and that's enough. I'm not in it for the big win at the end, whatever that means. I'm in it for what I find along the way.
 
 
Current Mood: optimisticoptimistic
Current Music: Wendy Carlos-Tron soundtrack