Computing the Fourier series of a matrix

It’s funny how you never know what kind of knowledge you are going to need in your work and life in general. For instance, the Spring I was making my Master’s Thesis (-13), I casually took a course on data structures and algorithms: stacks and queues and trees and the like. I held no great interest in it, and thought I would probably never need the information.

But lo’ and behold: the very second task I set my hands on the next Autumn as a reaaally junior Ph.D. student was the implementation of a reluctance network solver. And what else could that require than forming a number of fundamental loops for a graph by some spanning trees . In other words, exactly the stuff I thought I was never going to need.

Indeed, nowadays I try to stop myself from having thoughts like that.

…and I quite succesfully have. Nevertheless, I still regularly encounter stuff I never even knew existed. I mean algorithms totally unrelated to my field of studies, mathematical approaches almost forgotten, or simply a new search term to open up a world of new possibilities. Usually, these discoveries are related to solving some really minor problem of mine; just something that I need to get out of the way to move forward.

This time, I discovered highly oscillatory integrals. In other words, how to numerically integrate the product of a known function and a high-frequency sinusoidal function over a predetermined interval. This kind of integrals well be needed, for example, to calculate the fuckzillionth Fourier-coefficient of the aforementioned function.

And that is precisely what I’m doing, only with a highly unorthodox function. Indeed, I want to compute the Fourier-series of a matrix, specifically the moving-band matrix used in finite element (FE) simulations. And what would that be?

Well, when we are analysing a rotating electrical machine, we often mesh both the stator and rotor separately, and then connect these two meshes together with an air-gap mesh. Now, as the machine rotates, only this air-gap mesh has to be changed, yielding a significant speed-up for the computations. Typically, the mesh is regular enough that his movement can be handled by only changing one layer of the air-gap mesh: hence the name moving band.

In the end, the system to be solved can be written as a matrix equation
,
where is the vector of unknowns we are solving, and is the vector of known source-terms. Since we are dealing with a rotating system, both vectors will be time-dependent. Finally, the matrices and denote the static part of the system matrix, and the moving-band matrix related to the air-gap mesh, respectively.

Now, this moving band matrix is what I want to express as a Fourier-series. To do that, I need to calculate lots of integrals of type
,
one for each matrix entry and each angular frequency considered. The problem is made all the funnier by the facts that
1) There’s no simple analytical expression for the matrix : we can only obtain a snapshot of it at some instant of time .
2) The entries are only piecewise continuous in time, with discontinuities in between.
3) The matrix is somewhat computationally expensive to construct.

Due to 1) we need to calculate the Fourier-integral numerically, by sampling the matrix at different instants of time. This, in turn, is made somewhat difficult by 2) and 3). However, the presence of the oscillating term turns out to be even more problematic than this.

Indeed, below you can see some of the entries as functions of the rotor position. The blue one is some diagonal term, always non-zero. The others are some cross-coupling terms between the stator and rotor, only different from zero when they are close to each other. As can be seen, they are all far from smooth. Still, this situation could be handled by considering shorter sub-intervals on which the entries are both continuous and smooth, indicated by the vertical black dot-lines. On each sub-interval, the integral could then be calculated by e.g. quadratures.

…except this approach won’t really work with high-order terms, i.e. with . Due to the cosine function, the product terms start to look really ugly really soon. This is illustrated below for the relatively benign harmonic of order 49. As can be seen, the product terms start to get quite oscillatory indeed even in this case. In reality, terms up to the 1000th harmonic might be needed. As can be guessed, evaluating an integral of this kind numerically would require a really high number of sample points, resulting in nightmarish computation times.

Luckily, some fine fellow by the name of Filon had already come up with a solution for me in the early 20th century. His approach includes chopping the integration interval into smaller sub-intervals  – something I’d have to do in any case the handle the discontinuities. In each of these smaller intervals, the function (in my case ) would then be approximated with a second-order polynomial. Finally, the product of this polynomial and the cosine function could then be integrated in closed form without any problems with Wolfram Alpha.

I have to admit, this approach is so simplistically genious I’m ashamed I couldn’t figure it out myself. Well, you can’t win every time I guess. Anyway, I have already implemented it in Matlab, and it seems to work really well and surprisingly fast.

Nevertheless, I’m currently writing a C++ implementation of it, with the Eigen library. Since it’s been forever since I used any lower-level programming language, and have used C++ exactly never, things are not moving very fast right now. Still, I’m slowly inching forward, and hope that the resulting dll will be faster than my current Matlab implementation. Well, even if it isn’t I’ll still get some always-needed programming experience.

I’d call that a win.

Computing the Fourier series of a matrix