ApproxFun.jl Convolution/Fourier Multiplier

Specific Domains Numerics
1

I’m trying to implement a solver for the Benjamin-Ono equation using DifferentialEquations.jl. I have periodic boundary conditions, so I’m using ApproxFun.jl to take the Fourier transform and use the method of lines to solve the PDE.

The equation is defined as,

u_t + uu_x - Hu_{xx} = 0,

where H is the Hilbert transform defined as,

\hat{H} = -i \operatorname{sgn}(k)

I poked around the forums, GitHub examples, and the benchmark pages on DifferentialEquations.jl, but I wasn’t able to see an example of where someone used a linear operator other than a derivative. I’m a little confused about how ApproxFun.jl is storing the Fourier coefficients, and if I can get a grasp on that I shouldn’t have trouble implementing my own operator for the Hilbert transform.

Thanks!

EDIT:
Just to be a little more clear, I want the operator to be diagonal and represent the Fourier multiplier,

\mathcal{K} = i(k)^2 \operatorname{sgn}(k)

This way I bring the second derivative inside of the multiplier.

2

Fourier stores the cosine and sine coefficients alternating. Laurent stores the complex exponential coefficients alternating positive and negative coefficients, that is, 0, -1, 1, -2, 2, … coefficients.

The docs should explain in detail

2 Likes
3

Now looking at the docs this makes sense, I was unfamiliar with the notation that they use for the ordering. I’ve been working with the complex coefficients for so long it felt weird going back to the trigonometric coefficients :sweat_smile: Thanks for clearing up the notation in the docs for me.

Does anyone happen to know the benefit of working with ApproxFun.jl instead of the FFT implementation directly?

4

There’s not much benefit with fixed discretisation sizes, it’s more useful when you take advantage of the adaptivity. Or with non-periodic problems where operators are not diagonal.

I’d recommend using the trigonometric version: it ensures the solution remains real which may also be significantly faster. There is of course a simple expression for the Hilbert transform that swaps the sine and cosine coefficients.

1 Like
5

I have a follow-up question, I’m playing around with the trigonometric version, but I am also interested in the Taylor version. Glancing at the source code in the ApproxFun.jl package I was having trouble figuring out if it was using the plan_rfft, when playing around in Jupyter, it looked like the Taylor version was just using plan_fft. Is Taylor assuming that \hat{f}(k) = \hat{f}^*(-k), or is there something else going on?

6

A Taylor expansion is the same as Fourier with no negative terms. In this case, the “negative” coefficients of the FFT accurately approximate the large positive coefficients.

1 Like
7

Yes, prior to using DifferentialEquations.jl I was using an implementation of the symplectic Yoshida sixth-order integrator to solve PDEs. The adjoint PDE for Benjamin-Ono is non-atuonomous, so I’m trying to use some different solvers inside of the DifferentialEquations.jl package now. I’m not sure why, but I can’t seem to get the package to work with the FFT plans directly, so now I’m using ApproxFun.jl since this is what I’ve seen done in examples.

When I use the FFT I typically use the real FFT plan which uses conjugate symmetry to work with only the positive non-negative coefficients; like what the Taylor implementation is doing. It seems that it is not using the rfft_plan:

S = Taylor()
n = 512
T = ApproxFun.plan_transform(S, n)

typeof(T)

ApproxFunBase.TransformPlan{Float64,Hardy{true,Circle{Float64,Float64,Complex{Float64}},Complex{Float64}},false,ApproxFunBase.TransformPlan{Complex{Float64},Hardy{true,Circle{Float64,Float64,Complex{Float64}},Complex{Float64}},false,ApproxFunBase.TransformPlan{Complex{Float64},Hardy{true,Circle{Float64,Float64,Complex{Float64}},Complex{Float64}},true,FFTW.cFFTWPlan{Complex{Float64},-1,true,1}}}}

When I create a plan with the plan_rfft:

M = 512
u = Array{Float64}(undef, M)
T = plan_rfft(u)

typeof(T)

FFTW.rFFTWPlan{Float64,-1,false,1}

Is there a reason why Taylor is using cFFTWPlan instead of rFFTWPlan?

8

I don’t know why you’d want rfft… are you assuming extra symmetry?

In any case, rfft is not that much faster than fft. In fact it depends on the number of points, sometimes FFT is faster.

9

I don’t understand what additional symmetry is assumed?

The only assumption that rfft and Taylor is making is that the solution is real, which implies \hat{f}(k) = \hat{f}^*(-k). The r2r plan makes the additional assumption that the solution is also even, which implies \hat{f}(k) = \hat{f}(-k). However, if Taylor is also making conjugate symmetry assumption, then it doesn’t make sense that it’s not taking advantage of the rfftPlan. The reason this is preferable is because it requires half the space as the cfftPlan to store \hat{f}.

10

You can’t have real valued Taylor series… this is on the unit circle

11

Okay, this is what I was confused about. It is very possible to work with the series of just the positive Fourier coefficients, but you use the fact that \hat{f}(k) = \hat{f}^*(-k) for real solutions.

I work with real-periodic solutions. In these situations it is easier to only work with the non-negative Fourier coefficients. For me it is restrictive to only work with the trigonometric coefficients, and not have the luxury of taking advantage of symmetry properties of the complex coefficients. I’m struggling to understand the limitations of this package, and I appreciate your help.

12

That’s not a Taylor series. A Taylor series is an expansion f_0 + f_1 z + f_2 z^2 + .... There’s no support at the moment for Laurent series that takes advantage of symmetry in the coefficients.