100% Guaranteed Results


18-330 – Exercise 1: Fourier series Solved
$ 20.99
Category:

Description

5/5 – (1 vote)

In lectures we saw that we can write a periodic function ๐‘“ in the form
๐‘
๐‘“(๐‘ฅ) = โˆ‘ ๐‘“๐‘›ฬ‚ ๐‘’๐‘–๐‘›๐‘ฅ
๐‘›=โˆ’๐‘
where ๐‘“๐‘›ฬ‚ โˆถ= 21๐œ‹ โˆซ02๐œ‹ ๐‘“(๐‘ฅ)๐‘’โˆ’๐‘–๐‘›๐‘ฅd๐‘ฅ if the period is 2๐œ‹.
The FFTW library uses a different convention, taking the period to be 1. Suppose we sample a function at ๐‘ + 1 uniformly-spaced points between [0, 1], obtaining values ๐‘“0, ๐‘“1, โ€ฆ , ๐‘“๐‘, where ๐‘“๐‘› โˆถ= ๐‘“(๐‘›โ„Ž) with โ„Ž โˆถ= ๐‘1 . By using the Fourier series we are implicitly assuming that ๐‘“ is periodic, i.e. ๐‘“0 = ๐‘“๐‘. Then
๐‘2 โˆ’1
๐‘“(๐‘ฅ) = โˆ‘ ๐‘“๐‘›ฬ‚ ๐‘’2๐œ‹๐‘–๐‘›๐‘ฅ
๐‘›=โˆ’๐‘2
where the Fourier coefficients are
1
๐‘“๐‘›ฬ‚ โˆถ= โˆซ ๐‘“(๐‘ฅ) ๐‘’โˆ’2๐œ‹๐‘–๐‘›๐‘ฅd๐‘ฅ
0
This is the convention we will use in this problem.
1. Discretize the integral for the Fourier series coefficients using the trapezoidal rule using ๐‘ + 1 points. Using the assumption that ๐‘“0 = ๐‘“๐‘ reduce this to a sum over ๐‘“๐‘– for ๐‘– = 0, โ€ฆ , ๐‘ โˆ’ 1.
Implement this as a function fourier_coefficients_trapezoidal(f::Vector) where f is the vector of samples of ๐‘“๐‘› for ๐‘› = 0, โ€ฆ , ๐‘ โˆ’ 1. Your function should return the coefficients for ๐‘“๐‘˜ฬ‚ for ๐‘˜ = โˆ’๐‘2 , โ€ฆ , ๐‘2 โˆ’ 1 as a vector.
2. Consider the following three functions:
1. The sawtooth function
๐‘“(๐‘ฅ) = mod(๐‘ฅ, 1)
(The function mod(๐‘ฅ, 1) returns just the fractional part of a number.)
2. The โ€œWโ€ wave function
โŽง{๐‘ฅ 0 โ‰ค ๐‘ฅ < 0.5
๐‘”(๐‘ฅ) = 1 โˆ’ ๐‘ฅ 0.5 โ‰ค ๐‘ฅ < 1
โŽจ{โŽฉ๐‘”(mod(๐‘ฅ, 1)) otherwise
3. The function โ„Ž(๐‘ฅ) = exp(cos(2๐œ‹๐‘ฅ))
Plot them.
Given the properties of these functions, how would you expect their Fourier coefficients to decay?
3. Take ๐‘ = 100. Calculate ๐‘“,ฬ‚ ๐‘”ฬ‚and โ„Žฬ‚ using your fourier_ceofficients_trapezoidal function. Plot the magnitude of the coefficients as a function of ๐‘› only for ๐‘› โ‰ฅ 0. Do you see the expected behavior?
4. How accurately does our function fourier_coefficients_trapezoidal(f::Vector) calculate the Fourier coefficients? Use the following analytic solutions to calculate and plot the respective errors:
๐‘“0ฬ‚ = 12 and ๐‘“๐‘›ฬ‚ = 2 ๐œ‹๐‘›๐‘– for ๐‘› โ‰  0
1 1
๐‘”ฬ‚0 = 4 ; ๐‘”ฬ‚๐‘› = โˆ’๐œ‹ 2๐‘›2 for ๐‘› odd; and ๐‘”ฬ‚๐‘› = 0 for ๐‘› even
โ„Žฬ‚๐‘› = ๐ผ1(๐‘›)
where ๐ผ1(๐‘›) is the modified Bessel function of the first kind, calculated in Julia as besseli(n, 1) using the SpecialFunctions.jl package.
Calculate the error (using norm) as the number of points used changes between ๐‘ = 10 and ๐‘ = 1000.
5. Write a function reconstruct_fourier_series(fฬ‚s::Vector, xs::Vector) which reconstructs f(x) from the Fourier series coefficients.
Make a plot of ||f(xt) – reconstruct_fourier_series(fฬ‚s, xt)||/N as the number of coefficients used ๐‘ is varied from 4 โ†’ 200 for each of the functions.
Sample the reconstructed Fourier series at the points xt = 0:0.001:1. What do you see? Can you explain why?
6. Make an interactive visualiztion that plots the following on the same axes:
(i) the points used to calculate the Fourier coefficients in fourier_coefficients_rectangle;
(ii) the reconstructed function from the Fourier coefficients, found using reconstruct_fourier_series; and
(iii) the true function as the number of points is varied from ๐‘ = 10 โˆถ 2 โˆถ 250.
Does this help explain the results in [1.6]? In particular, what do you see for the sawtooth function? This is known as the Gibbs phenomenon, which occurs when a function is discontinuous.
7. What is the operation count for your naive fourier_coefficients_trapezoidal function? In general this will not behave well as ๐‘ grows very large. The FFT, however, can calculate this in ๐‘‚(log(๐‘)) steps. The FFT implemented in FFTW.jl calculates
๐‘โˆ’1
๐‘“๐‘›ฬ‚ = โˆ‘ ๐‘“๐‘˜๐‘’โˆ’2๐œ‹๐‘–๐‘›๐‘˜/๐‘
๐‘˜=0
for ๐‘› = 0 โˆถ ๐‘ โˆ’ 1, which should be related to your discretization in [1.1].
Note that the ๐‘›s are different here. But since ๐‘’โˆ’2๐œ‹๐‘–(๐‘›+๐‘™๐‘)๐‘˜/๐‘ =
๐‘’โˆ’2๐œ‹๐‘–๐‘™๐‘˜๐‘’โˆ’2๐œ‹๐‘–๐‘›๐‘˜/๐‘ = ๐‘’โˆ’2๐œ‹๐‘–๐‘›๐‘˜/๐‘ we have the relationship ๐‘“๐‘/2+๐‘–ฬ‚ =
๐‘“โˆ’ฬ‚ ๐‘/2+๐‘–.
The FFT algorithm therefore outputs
fฬ‚= [๐‘“0ฬ‚ , ๐‘“1ฬ‚ , โ€ฆ , ๐‘“๐‘/2โˆ’1ฬ‚ , ๐‘“โˆ’ฬ‚ ๐‘/2, ๐‘“โˆ’ฬ‚ ๐‘/2+1, โ€ฆ , ๐‘“โˆ’1ฬ‚ ]
The FFTWpackage defines a function fftshift that shifts this vector to the form
fฬ‚= [๐‘“โˆ’ฬ‚ ๐‘/2, ๐‘“โˆ’ฬ‚ ๐‘/2+1, โ€ฆ , ๐‘“โˆ’1ฬ‚ , ๐‘“0ฬ‚ , ๐‘“1ฬ‚ , โ€ฆ , ๐‘“๐‘/2โˆ’1ฬ‚ ]
8. Implement a function fast_fourier_coefficients that outputs the same results as fourier_coefficients_trapezoidal but using the FFT from FFTW.jl.
Check that the output is the similar to before.
Time your two functions for ๐‘ = 210. Is one faster than the other?
How large can you take N such that fast_fourier_coefficients runs for under 1 second? What about for fourier_coefficients_trapezoidal?
Exercise 2: Solving an ODE with a spectral method
In this problem we will solve the boundary-value problem
๐‘“โ€ณ = ๐‘”
with boundary conditions ๐‘“(0) = ๐‘“(1) = 0 by using a spectral method, i.e. by expanding in suitable basis functions.
๐‘“(๐‘ฅ) = (๐‘›๐œ‹๐‘ฅ),
which is the same as the Fourier series for an odd function.
In practice we have to truncate the summation as
๐‘
๐‘“(๐‘ฅ) = โˆ‘ ๐‘๐‘› sin(๐‘›๐œ‹๐‘ฅ)
๐‘›=1
The coefficients are given by
1
๐‘๐‘› = 2 โˆซ ๐‘“(๐‘ฅ) sin(๐‘›๐œ‹๐‘ฅ)d๐‘ฅ
0
Similarly to the DFT we discretize this using the rectangle rule to get the Discrete Sine Transform (DST): ๐‘โˆ’1
2 ๐‘›๐‘˜๐œ‹
๐‘๐‘› = ๐‘“(๐‘ฅ) sin ( ๐‘ )
This sum is implemented in julia using the FFTW.jl library (which you will need to install) as follows:
dst(x) = FFTW.r2r(x, FFTW.RODFT00) / length(x) Its inverse is given by idst(x) = FFTW.r2r(x, FFTW.RODFT00)
[r2r stands for โ€œreal to realโ€, meaning that the transform maps a real vector to a real vector. RODFT00 is a symbol that selects one particular type of transform.]
1. Assume that there is a solution of the form ๐‘“(๐‘ฅ) = โˆ‘ฬƒ ๐‘๐‘›=1 ๐‘๐‘› sin(๐‘›๐œ‹๐‘ฅ).
Substitute this into the ODE ๐‘“โ€ณ(๐‘ฅ) = ๐‘”(๐‘ฅ) to show that ๐‘“๐‘›ฬ‚ฬƒ ๐œ‹ ๐‘› .
2. We can therefore solve the ODE for ๐‘“ by first calculating ๐‘”ฬ‚using the DST, then calculate ๐‘“๐‘›ฬ‚ , and finally invert using the iDST.
Write a function spectral_solver(b) that does this to solve the ODE, where b is the discretized version of g(x). Solve the ODE with ๐‘”(๐‘ฅ) = sin(2๐œ‹๐‘ฅ). Plot the result. The right-hand side is given by
h = 1 / N
b = sin.(2ฯ€ * (h:h:1-h))
3. Calculate the error as a function of ๐‘ and plot it. What rate of convergence do you see?
4. Generate the error plots again, now for the right-hand side ๐‘”(๐‘ฅ) = exp(sin(2๐œ‹๐‘ฅ)) โˆ’ 1. Use the solution from your spectral solver with
๐‘ = 213 as the true solution. Calculate the error for ๐‘ = 2๐‘› for ๐‘› = 1 โ†’ 12. Some care is needed to make sure you use the correct points from the โ€œtrueโ€ solution for comparison.
Exercise 3: Finding roots in a different way
In class we defined the Chebyshev expansion of a function as
๐‘
๐‘“๐‘(๐‘ฅ) = โˆ‘ ๐‘๐‘›๐‘‡๐‘›(๐‘ฅ)
๐‘›=0
which is an ๐‘th-degree polynomial. The Chebyshev polynomials are defined as ๐‘‡๐‘›(๐‘ฅ) = cos(๐‘› arccos(๐‘ฅ)).
In general for a smooth function the Chebyshev series converges rapidly. We therefore expect that the roots of ๐‘“๐‘(๐‘ฅ) should be close to the roots of ๐‘“(๐‘ฅ), provided that ๐‘“๐‘ is indeed a good approximation to ๐‘“. We have already seen that we can find all the roots of a polynomial with various methods.
The Chebyshev polynomials satisfy the following recurrence relation:
๐‘‡๐‘›+1(๐‘ฅ) = 2๐‘ฅ ๐‘‡๐‘›(๐‘ฅ) โˆ’ ๐‘‡๐‘›โˆ’1(๐‘ฅ)
We will use this to find the companion matrix for ๐‘“๐‘(๐‘ฅ), from which we can find the roots of ๐‘“๐‘(๐‘ฅ).
1. Consider the polynomial ๐‘ฅ๐‘‡๐‘›(๐‘ฅ). This is a degree ๐‘›+1 polynomial and hence can be reexpanded in Chebyshev polynomials.
Consider the vector of Chebyshev polynomials
๐‘‡0(๐‘ฅ)
T(๐‘ฅ) = โŽข 1(โ‹ฎ๐‘ฅ) โŽคโŽฅ โŽก ๐‘‡
โŽฃ๐‘‡๐‘โˆ’1(๐‘ฅ)โŽฆ
Now we can write ๐‘ฅT(๐‘ฅ) = ๐‘€T(๐‘ฅ) + ๐‘˜๐‘๐‘’ฬ‚๐‘๐‘‡๐‘(๐‘ฅ), where ๐‘€ is an (๐‘ โˆ’1)ร—(๐‘ โˆ’1) matrix. We need the ๐‘˜๐‘ term to account for the fact that ๐‘ฅ๐‘‡๐‘โˆ’1 is a degree-๐‘ polynomial. Here, ๐‘’ฬ‚๐‘ is the standard basis vector with zeros everywhere except in the ๐‘th component.
Use the recurrence relation to find the form of ๐‘€ and ๐‘˜๐‘. 2. Verify what you found in [3.1] numerically when ๐‘ = 10. Build the vector T(๐‘ฅ) for ๐‘ฅ a random number in [โˆ’1, 1]. Compute ๐‘€ and check that it gives ๐‘ฅT(๐‘ฅ).
3. This looks almost like an eigenvalue problem, except for the ๐‘‡๐‘ term. To remove this we can add and subtract ๐ถ๐‘“๐‘(๐‘ฅ)๐‘’ฬ‚๐‘ from the right-hand
side. Writing ๐‘“๐‘(๐‘ฅ) = ๐‘๐‘๐‘‡๐‘(๐‘ฅ) + โˆ‘๐‘โˆ’1๐‘›=0 ๐‘๐‘›๐‘‡๐‘›(๐‘ฅ), what value of ๐ถ should you choose so that
๐‘โˆ’1
๐‘ฅT(๐‘ฅ) = ๐‘€T(๐‘ฅ) + ๐ถ๐‘“๐‘(๐‘ฅ)๐‘’ฬ‚๐‘ โˆ’ ๐ถ๐‘’ฬ‚๐‘ โˆ‘ ๐‘๐‘›๐‘‡๐‘›(๐‘ฅ) (1)
๐‘›=0
= ๐‘€ฬƒ T(๐‘ฅ) + ๐ถ๐‘“๐‘(๐‘ฅ)๐‘’ฬ‚๐‘ (2)
What is the new matrix ๐‘€ฬƒ ?
4. This becomes an eigenvalue problem when ๐‘ฅ is a root of ๐‘“๐‘(๐‘ฅ). Therefore, the eigenvalues of ๐‘€ฬƒ are the roots of ๐‘“๐‘.
Write a function buildM(c::Vector) that constructs the matrix ๐‘€ฬƒ from the coefficients in the Chebyshev expansion. Use this to write a chebyshev_roots(c) function that finds the roots of the polynomial defined using the Chebyshev coefficients c. Finally write a function fN(x, c) that calculates the series expansion to find ๐‘“๐‘(๐‘ฅ) defined by the vector c. 5. We can calculate Chebyshev coefficients using the dct functions in FFTW.
We will use the Chebyshev points ๐‘ฅ๐‘› โˆถ= cos(๐‘›๐œ‹/๐‘) for ๐‘› = 0 โˆถ ๐‘. You can then calculate the Chebyshev coefficients using the following code:
chebyshev_points(N) = cos.(ฯ€*(0:1:N)/N)
function chebyshev_coefficients(x) N = length(x)
c = FFTW.r2r(x, FFTW.REDFT00)/(N-1)
c[1] /= 2 c[N] /=2
return c
end 6. Consider the polynomial ๐‘“(๐‘ฅ) = ๐‘ฅ (๐‘ฅ โˆ’ 1/2)2 (๐‘ฅ2 โˆ’ 1/9) (๐‘ฅ + 1/4).
Using 10 Chebyshev points calculate the Chebyshev coefficients and then calculate the roots using chebyshev_roots. What do you see. What about multiplicities?
7. Plot ๐‘“(๐‘ฅ) and scatter plot the roots you find on top.
8. Now consider solving the problem exp(cos(4๐œ‹๐‘ฅ)) = 1. Using ๐‘ = 100 points, alculate the Chebyshev coefficients for ๐‘”(๐‘ฅ) = exp(cos(4๐œ‹๐‘ฅ))โˆ’ 1. Do they decay quickly? Use these to calculate the roots of ๐‘”(๐‘ฅ). Plot ๐‘”(๐‘ฅ) and scatter the roots you find on top. Do you find all the roots?
(Hint: you will find 100 eigenvalues. Only plot those that are real and lie between -1 and 1.)
9. Make an interactive visualization as ๐‘ is varied between ๐‘ = 4 โˆถ 150. Plot ๐‘”(๐‘ฅ), the Cheyshev approximation to ๐‘”(๐‘ฅ) using ๐‘ coefficients and the roots you find on the same axes. Comment on what you see.
At what value of ๐‘ do you find all the roots? Remember to plot only those roots that are real and between -1 and 1.
Exercise 4: Gramโ€“Schmidt for polynomials In lectures we discussed treating the set of polynomials {1, ๐‘ฅ, ๐‘ฅ2, ๐‘ฅ3, โ€ฆ} as the basis of a vector space with the inner product
๐‘
(๐‘“, ๐‘”)๐‘ค = โˆซ ๐‘“(๐‘ฅ)๐‘”(๐‘ฅ)๐‘ค(๐‘ฅ)d๐‘ฅ
๐‘Ž
We can therefore carry out Gramโ€“Schmidt orthogonalization on these polynomials to generate a family of orthogonal polynomials.
We will implemnent this using the Polynomials.jl package. Integrals can be performed using the polyint(f, a, b) function to integrate the polynomial ๐‘“ from ๐‘Ž to ๐‘. Here, ๐‘“, ๐‘” and ๐‘ค should all be Polynomials.
1. Write a function gram_schmidt(vs::Vector, ip) which accepts a vector of โ€œvectorsโ€ in the vector space and the inner product on the vector spacce. For standard vectors this would be dot(v1, v2). The function should implement the Gramโ€“Schmidt algorithm and return a vector of the resulting orthonormal basis elements.
2. Test your function for standard vectors vs = [rand(10) for i = 1:10] using ip = dot. To check everything went according to plan, form the matrix ๐ท๐‘๐‘–๐‘— = ๐‘ž๐‘– โ‹… ๐‘ž๐‘—. If everything worked this should be the identity matrix.
3. For polynomials we define the inner product
๐‘
(๐‘“, ๐‘”)๐‘ค โˆถ= โˆซ ๐‘“(๐‘ฅ)๐‘”(๐‘ฅ)๐‘ค(๐‘ฅ)d๐‘ฅ
๐‘Ž
For Legendre polynomials we have ๐‘Ž = โˆ’1 , ๐‘ = 1, and ๐‘ค(๐‘ฅ) =
1. Using the Polynomials.jl package, implement a function legendre_inner_product(f,g). Use this to orthogonalize the vector of monomials up to order 7.
4. Use the functions you found in [4.3] and your inner-product function, find the Legendre polynomial expansion coefficients for the function ๐‘“(๐‘ฅ) = ๐‘ฅ (๐‘ฅ โˆ’ 1/2)2 (๐‘ฅ2 โˆ’ 1/9) (๐‘ฅ + 1/4).
Plot the reconstructed polynomial using the first ๐‘› coefficients for ๐‘ = 1 โˆถ 7 and the true function. What do you see? How good are the Legendre polynomials at approximating the function?

Reviews

There are no reviews yet.

Be the first to review “18-330 – Exercise 1: Fourier series Solved”

Your email address will not be published. Required fields are marked *

Related products