@@ -15,11 +15,11 @@ n = 256
1515L = 2π
1616
1717dx = L / n
18- x = range(start=-L/ 2, stop=L/2- dx, length= n) |> Array
19- u = @. sin(5x)cos(7x);
18+ x = range(start = -L / 2, stop = L / 2 - dx, length = n) |> Array
19+ u = @. sin(5x)cos(7x);
2020du = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
2121
22- k = rfftfreq(n, 2π*n/ L) |> Array
22+ k = rfftfreq(n, 2π * n / L) |> Array
2323m = length(k)
2424P = plan_rfft(x)
2525
@@ -29,23 +29,19 @@ bwd(u, p, t) = P \ u
2929fwd(du, u, p, t) = mul!(du, P, u)
3030bwd(du, u, p, t) = ldiv!(du, P, u)
3131
32- F = FunctionOperator(fwd, x, im*k;
33- T=ComplexF64,
34-
35- op_adjoint = bwd,
36- op_inverse = bwd,
37- op_adjoint_inverse = fwd,
38-
39- islinear=true,
40- )
32+ F = FunctionOperator(fwd, x, im * k;
33+ T = ComplexF64, op_adjoint = bwd,
34+ op_inverse = bwd,
35+ op_adjoint_inverse = fwd, islinear = true
36+ )
4137
4238ik = im * DiagonalOperator(k)
4339Dx = F \ ik * F
4440
4541Dx = cache_operator(Dx, x)
4642
47- @show ≈(Dx * u, du; atol= 1e-8)
48- @show ≈(mul!(copy(u), Dx, u), du; atol= 1e-8)
43+ @show ≈(Dx * u, du; atol = 1e-8)
44+ @show ≈(mul!(copy(u), Dx, u), du; atol = 1e-8)
4945```
5046
5147## Explanation
@@ -60,14 +56,13 @@ FFT wrapper.
6056using SciMLOperators
6157using LinearAlgebra, FFTW
6258
63- L = 2π
64- n = 256
59+ L = 2π
60+ n = 256
6561dx = L / n
66- x = range(start=-L/ 2, stop=L/2- dx, length= n) |> Array
62+ x = range(start = -L / 2, stop = L / 2 - dx, length = n) |> Array
6763
68- u = @. sin(5x)cos(7x);
64+ u = @. sin(5x)cos(7x);
6965du = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x);
70-
7166```
7267
7368Now, we define the Fourier transform. Since our input is purely Real, we use the real
@@ -77,8 +72,8 @@ and `LinearAlgebra.mul!(xhat, transform, x)`. We also get `k`, the frequency mo
7772our finite grid, via the function ` rfftfreq ` .
7873
7974``` @example fft_explanation
80- k = rfftfreq(n, 2π*n/ L) |> Array
81- m = length(k)
75+ k = rfftfreq(n, 2π * n / L) |> Array
76+ m = length(k)
8277P = plan_rfft(x)
8378```
8479
@@ -93,15 +88,11 @@ bwd(u, p, t) = P \ u
9388
9489fwd(du, u, p, t) = mul!(du, P, u)
9590bwd(du, u, p, t) = ldiv!(du, P, u)
96- F = FunctionOperator(fwd, x, im*k;
97- T=ComplexF64,
98-
99- op_adjoint = bwd,
100- op_inverse = bwd,
101- op_adjoint_inverse = fwd,
102-
103- islinear=true,
104- )
91+ F = FunctionOperator(fwd, x, im * k;
92+ T = ComplexF64, op_adjoint = bwd,
93+ op_inverse = bwd,
94+ op_adjoint_inverse = fwd, islinear = true
95+ )
10596```
10697
10798After wrapping the FFT with ` FunctionOperator ` , we are ready to compose it with other
@@ -115,6 +106,6 @@ Dx = F \ ik * F
115106
116107Dx = cache_operator(Dx, x)
117108
118- @show ≈(Dx * u, du; atol= 1e-8)
119- @show ≈(mul!(copy(u), Dx, u), du; atol= 1e-8)
109+ @show ≈(Dx * u, du; atol = 1e-8)
110+ @show ≈(mul!(copy(u), Dx, u), du; atol = 1e-8)
120111```
0 commit comments