Improving the convergence of the Fourier-series expansion of ring currents for simulating nonaxisymmetric dipole sources in cylindrical coordinates

[oskooi]

The simulation of a nonaxisymmetric dipole source in cylindrical coordinates involves a Fourier-series expansion of "ring" current sources with azimuthal dependence $e^{im\phi}$ for $m = 0, 1, 2, ..., M$ for some cutoff $M$ which is proportional to the dipole position $r$. When the dipole is positioned far from $r = 0$, $M$ can become large (i.e., 50+) which can significantly increase the total runtime.

It would be useful to somehow speed up the convergence of the Fourier-series expansion.

[stevengj]

What you're looking for is some kind of series extrapolation method. Not sure which one works best for this particular problem, but you could try a few and experiment.

Basically what you need is a test problem where you've calculated it for lot's of $m$'s, and then compare it to plugging in the first few into a series-extrapolation code, and see which method converges the fastest.

[oskooi]

To experiment with the series extrapolation, here is the Fourier-series data for the source and radiated power of a ring current vs. $m$ for a nonaxisymmetric dipole (at $r$ = 6.3 μm) in a dielectric slab ($n$ = 2.4, thickness [in vacuum] of 0.7 um) above a perfect-reflector ground plane. This setup is based on the tutorial example. The script used to generate these results is here.

In this example, the cutoff $M$ is determined by when the radiated power (in air) drops to below 1% its maximum value over the interval of all preceding $m$ values. The value for $M$ found numerically is consistent with the analytic value of roughly $\omega r$ (~42).

The series data is provided at resolutions 50 and 100:

resolution_50.csv

resolution_100.csv

[stevengj]

You'll need a lot more m's before we can accurately brute-force the $m \to \infty$ limit (for comparison purposes with the extrapolations). Because the series is oscillatory, it's not enough to stop when a single $m$ goes below a 1% threshold — you neeed to look at the peak of the oscillations.

[stevengj]

(I'm not sure how well Richardson will work on extrapolating an oscillatory series like this — the oscillatory part seems like it will make a polynomial fit less accurate. I found a 2005 paper on a related problem that might have some suggestions, but it looks fairly empirical for one specific case that is different from ours.)

The Python mpmath.nsum function includes several extrapolation methods that would be worth looking into. Its default is a combination of Richardson and Shanks transformations, which might help here.

[stevengj]

There is a 2003 book by Sidi on extrapolation methods that may be helpful too.

Another thing to keep in mind is that we lose some information when we take the real part of the power — for extrapolation purposes, it might be better to keep both the real and imaginary parts, extrapolate, and take the real part at the end.

[stevengj]

I did some experiments with the slowly converging oscillatory series (derived from the Fourier series of a square window function), in fact only conditionally converging:

$\sum_{m=1}^{\infty} \frac{\sin m}{m} = \frac{\pi - 1}{2} \approx 1.0707963267948966$

whose partial sums

$S_n = \sum_{m=1}^{n} \frac{\sin m}{m}$

look like

Not surprisingly, Richardson extrapolation of $S_n$, treated as a function of $h = 1/n$ and extrapolating $h \to 0^+$, is terrible — it actually makes the error worse. That's because the oscillations cause the $S_n$ to look nothing like a polynomial in $1/n$. A Shanks transformation of $S_n$ also doesn't seem to help.

However, if I instead consider the complex partial sums:

$C_n = \sum_{m=1}^{n} \frac{e^{im}}{m}$

for which $S_n = \text{Im}(C_n)$, then a Shanks transformation helps a lot. (That is, you do a Shanks transformation of $C_n$ and then take the imaginary part.) I think the basic reason for this is that the Shanks transformation is based on an exponential ansatz, and the $e^{im}$ terms in $C_n$ are an exponential while the $\sin m$ terms in $S_n$ are not. Here is a plot of the convergence errors of the respective sums:

(Even after the Shanks transformation of $C_n$, Richardson extrapolation is still not helpful in this case.)

If this is at all similar to the types of Fourier series we are summing in Meep, keeping the full complex quantity and doing a Shanks transformation might be worthwhile.

[smartalecH]

@stevengj would a Pade approximate over the partial sums work?

[stevengj]

Pade might be another option, but again I think we would ideally want the complex coefficients.

I'm not certain of the best way to use Pade here. The most obvious way to use Pade would be to apply it the fields, which form a Fourier-series polynomial in $e^{i\phi}$, and use Pade to extrapolate them to the total fields. I'm less sure at first glance how to apply Pade extrapolation to something like the total power.