AMATH242/CS371 Spring 2024 Assignment IV
Release date: Tuesday, July 16th
Due date: Tuesday, July 30th, 11:59 pm
One-time extension not applicable
• Questions below are either theoretical or computational. For the computational questions you may use any pro- gramming language you prefer except symbolic ones like Maple or Mathematica.
• This assignment should be submitted to Crowdmark. Besides your written/typed solutions you must also submit your code. Please make sure your upload has the correct orientation and ordering.
• If you want to use your one time 3-day extension, please email the instructor before the original due date.
• You are not allowed to post this assignment on sites like stackexchange.com, chegg.com, etc. This will be checked.
Offenders will face penalties from the Math Faculty (e.g. suspension).
Total points: 20
1. (0points)Please sign the Academic Integrity Checklist on the lastpage ofthis pdf. Ifyou do not sign the Academic Integrity Checklist you will receive a 0 for this assignment.
2. (Theoretical, 3 points) Given {(x0 , f0 ), (x1 , f1 ), (x2 , f2 ) }, we want to find a degree 4 spline on [x0 , x2 ]. To deter- mine the spline, we need three groups of conditions: usual interpolation conditions; smoothness conditions and extra conditions. You may use p[1] (x) to denote the part of the spline on [x0 , x1 ] and p[2] (x) to denote the part of the spline on [x1 , x2 ].
(1) How many unknowns do we have in this spline? For example, if we want to determine a polynomial a0 + a1 x + a2 x2 , there are, in this case, three unkowns a0 , a1 , a2 .
(2) List out all usual interpolation conditions and smoothness conditions.
3. (Computational, 8 points) We seek to numerically compute the double integral:
∫ ∫ sin(x) cos(y) dxdy. (♣)
The quadrature rules covered in [YW] integrate a function of one variable between two points. To compute eq. (♣), we observe that
sin(x) cos(y) dxdy = sin(x) dx) . cos(y) dy) .
We deduce that the exact value of eq. (♣) is (1 − cos(1))sin(1). Also, we can construct a quadrature rule for eq. (♣) by computing the following:
whereI(̂) can be any quadrature rule that is covered in [YW].
We use composite quadrature rules in this question. As usual, we divide [0, 1] into equal subintervals, with each subinterval’s length denoted by h. This applies to both sin(x) dx and cos(y) dy.
(1) Implement composite trapezoidal rule for eq. (♣). In other words, computeI(^)CT(sin(x)) . I(^)CT(cos(x)). Your
program is required to do the following:
• Compute I(^)CT(sin(x)) . I(^)CT(cos(x)) for h = 1/5, 1/10, 1/20, 1/40, 1/80.
• For each h, compute the error:
E = (1 − cos(1)) sin(1) − I(^)CT(sin(x)) . I(^)CT(cos(x)).
• Compute the order of convergence. To do this, take two consecutive h’s, say hi−1, hi, and their corre- sponding errors Ei−1,Ei and compute:
Compute pi for all consecutive pairs of (h,E)’s. Note that when the error behaves like E = o(hp), these computed pi’s should be close to p.
Tabulate your results as follows:
|
i
|
1
|
2
|
3
|
4
|
5
|
|
hi
|
1/5
|
1/10
|
1/20
|
1/40
|
1/80
|
|
|Ei|
|
|
|
|
|
|
|
pi
|
N/A
|
|
|
|
|
(2) Implement composite Simpson’s rule to approximate eq. (♣). In other words, compute I(^)CS(sin(x)).I(^)CS(cos(x)).
Repeat all the steps required in the previous subquestion and tabulate your final results in the same manner.
(3) In [YW], we covered that|ECT | ≤ o(h2) and|ECS | ≤ o(h4). Show that if|ECT | = o(h2) and|ECS | = o(h4), then
You may assume that f(x) dx, g(y) dy,I(^)CT(f),I(^)CS(f),I(^)CT(g),I(^)CS(g) are all o(1) quantities (i.e., some
nonzero constants independent of h).
Make a brief 1-2 sentence comment relating the theoretical estimates with your computed pi’s from sub- questions (1) and (2).
4. (Theoretical, 4 points) We start with a real-valued signal vector of length 2N:
{f[n]}n(N)=−N+1 = {f[−N + 1],f[−N + 2], … ,f[−1],f[0],f[1], … ,f[N − 1],f[N]} .
Like what we did in the lecture notes, we assume that the signal has period 2N, which implies that f[−N] = f[N]. We further assume that the signal has odd symmetry, i.e. f[n] = −f[−n].
(1) Show that f[0] = f[N] = 0.
(2) Apply the DFT to {f[n]}n(N)=−N+1 and simplify the resulting formula for F[k] to the following form.
for k = −N + 1 ∼ N.
(3) Are F[k]’s all real-valued or all purely imaginary, or neither is true?
(4) Are {F[k]}k(N)=−N+1 even in k, or odd in k or neither is true?
5. (Computational, 5 points) Figure1shows a plot of a noisy signal (in blue). The signal shows both a pattern and irregularities. The pattern - a smooth signal - is drawn as a red curve while irregularities (oscillations) are due to added Gaussian noise.
The data file noisy signal.txt on LEARN provides the noisy signal (in blue) vector. To import the data and plot the signal in MATLAB, you may do:
Figure 1: A signal which is polluted by Gaussian noise.
1 x = load (' noisy_signal .txt ') ;
2 plot ( x ) ;
which generates the noisy signal (in blue) as in fig. 1. In Python, data can be imported and the same plot can be genrated with the following script.
1 import numpy as np
2 import matplotlib . pyplot as plt
3 x = np . loadtxt (' noisy_signal .txt ')
4 plt . plot ( x )
5 plt . show ()
(1) Implement and find the DFT of the noisy signal imported from noisy signal.txt using the following
|
formula:
|
|
|
|
|
0 ≤ k ≤ N − 1.
|
(♥1)
|
The output is a generally complex-valued vector {F [k]}k(N)01 . In MATLAB,keep in mind that the index needs
to be adjusted since MATLAB’s index starts from 1. For each k, compute the power which is given by:
N |F [k]|2 .
This quantity measures how strong a presence a certain mode k has in the signal. Make a plot with power as the y-axis and mode k as the x-axis. This plot is known as the power spectrum. Make brief observations on your power spectrum plot.
(2) Perform. the following tasks:
• Filter out the noise from the signal;
• Recover the smooth signal (the red curve in fig. 1);
• Make a plot of the denoised signal which should be close to the red curve in fig. 1.
For the first task, by trial and error, choose a suitable threshold of power and set every mode to zero when
it has power below the threshold, i.e.,
Set F [k] to 0, if N |F [k]|2 < .
For the second task, implement the IDFT using:
k , 0 ≤ n ≤ N − 1.
k=0
A useful sanity check of your code: apply eq. (♥2) to the output vector of eq. (♥1) without any filtering. Since eq. (♥2) is the inverse of eq. (♥1), you should recover the original signal vector. (You are not required to perform. this check in your submission.)
Some pointers in using MATLAB/NumPy with complex numbers:
• The imaginary unit is 1j in both MATLAB and NumPy;
• In NumPy, to initiate azero vector in complex number data type, you might need np.zeros(n,dtype=complex);
• After IDFT, you would obtain a complex-valued vector with imaginary part of each element being close to zero. Let’s say the vector’s name is f filtered, you might need to use plot(real(f filtered)) in
MATLAB and plt.plot(f filtered.real) in NumPy for correct plotting.