Fourier Analysis & Non-Linear Oscillations in Computational Physics: Lecture Notes - Prof., Study notes of Physics

Download Fourier Analysis & Non-Linear Oscillations in Computational Physics: Lecture Notes - Prof. and more Study notes Physics in PDF only on Docsity! PHYS-4007/5007: Computational Physics Course Lecture Notes Section XI Dr. Donald G. Luttermoser East Tennessee State University Version 4.1 Abstract These class notes are designed for use of the instructor and students of the course PHYS-4007/5007: Computational Physics taught by Dr. Donald Luttermoser at East Tennessee State University. Donald G. Luttermoser, ETSU XI–3 riodic =⇒ any function that repeats itself is said to be periodic. d) More precisely, if there exists a finite number L such that f(x + L) = f(x), then f(x) is periodic with period L. e) We can write any function that is periodic (or that is defined on a finite interval) as a Fourier series. f) However if f(x) is non-periodic or is defined on the infinite interval from −∞ to +∞, we must use a Fourier integral. 3. Fourier Series. Fourier series are not mere mathematical de- vices; they can be generated in the laboratory (or telescope) =⇒ a spectrometer decomposes an electromagnetic wave into spectral lines, each with a different frequency and amplitude (intensity). Thus, a spectrometer decomposes a periodic function in a fashion analogous to the Fourier series. a) Suppose we want to write a periodic, piecewise continuous function f(x) as a series of simple functions. Let L denote the period of f(x), and choose as the origin of coordinates the midpoint of the interval defined by this period−L/2 ≤ x ≤ L/2. b) If we let an and bn denote (real) expansion coefficients, we can write the Fourier series of this function as f(x) = a0 + ∞∑ n=1 [ an cos ( 2πn x L ) + bn sin ( 2πn x L )] . (XI-3) c) We calculate the coefficients in Eq. (XI-3) from the func- tion f(x) as a0 = 1 L ∫ L/2 −L/2 f(x) dx (XI-4) XI–4 PHYS-4007/5007: Computational Physics an = 2 L ∫ L/2 −L/2 f(x) cos ( 2πn x L ) dx (n = 1, 2, . . .) (XI-5) bn = 2 L ∫ L/2 −L/2 f(x) sin ( 2πn x L ) dx (n = 1, 2, . . .) (XI-6) d) Notice that the summation in Eq. (XI-3) contains an in- finite number of terms. In practice we retain only a finite number of terms =⇒ this approximation is called trun- cation. i) Truncation is viable only if the sum converges to whatever accuracy we want before we chop it off. ii) Truncation is not as extreme an act as it may seem. If f(x) is normalizable, then the expan- sion coefficients in Eq. (XI-3) decrease in magni- tude with increasing n, i.e., |an| → 0 and |bn| → 0 as n → ∞. iii) Under these conditions, which are satisfied by physically admissible wave functions, the sum in Eq. (XI-3) can be truncated at some finite maxi- mum value nmax of the index n. (Trial and error is typically needed to determine the value of nmax that is required for the desired accuracy.) iv) If f(x) is particularly simple, all but a small, fi- nite number of coefficients may be zero. One should always check for zero coefficients first before evalu- ating the integrals in Eqs. (XI-4, 5, 6). Donald G. Luttermoser, ETSU XI–5 4. The Power of Parity. One should pay attention as to whether one is integrating an odd or an even function. Trigonometric functions have the well-known parity properties: sin(−x) = − sin x (odd) (XI-7) cos(−x) = + cos x (even) (XI-8) As such, if f(x) is even or odd, then half of the expansion coef- ficients in its Fourier series are zero. a) If f(x) is odd [f(−x) = −f(x)], then    an = 0 (n = 1, 2, . . .) f(x) = ∞∑ n=1 bn sin ( 2πn x L ) (XI-9) b) If f(x) is even [f(−x) = +f(x)], then    bn = 0 (n = 1, 2, . . .) f(x) = ∞∑ n=1 an cos ( 2πn x L ) (XI-10) c) If f(x) is either an even or and odd function, it is then said to have definite parity. 5. The Complex Fourier Series: If f(x) does not have a definite parity, we can expand it in a complex Fourier series. a) To derive this variant on the Fourier series in Eq. (XI- 3), we just combine the coefficients an and bn so as to introduce the complex exponential function ei2πnx/L; viz., f(x) = ∞∑ n=−∞ cne i2πnx/L. (XI-11) b) Note carefully that in the complex Fourier series in Eq. (XI-11) the summation runs from −∞ to ∞. The expan- XI–8 PHYS-4007/5007: Computational Physics rewritten as H(ω) = ∫ ∞ −∞ h(t) eiωt dt (XI-19) h(t) = 1 2π ∫ ∞ −∞ H(ω) e−iωt dω . (XI-20) a) To introduce symmetry between these two equations, of- ten the 1/2π coefficient is split between the two integrals, introducing a 1/ √ 2π coefficient to each equation (as done in Eqs. XI-1,2,13,14). b) For this section, we will follow the f notation. 5. As we have seen, the following statements about these functions can be made: If ... then ... h(t) is real H(−f) = H∗(f) h(t) is imaginary H(−f) = −H∗(f) h(t) is even H(−f) = H(f) [i.e., even] h(t) is odd H(−f) = −H(f) [i.e., odd] h(t) is real and even H(f) is real and even h(t) is real and odd H(f) is imaginary and odd h(t) is imaginary and even H(f) is imaginary and even h(t) is imaginary and odd H(f) is real and odd These symmetries will be useful in order to develop computa- tional efficiency in coding. 6. Useful scalings and shifting equations: h(t) ⇐⇒ H(f) “no scaling” (XI-21) h(at) ⇐⇒ 1 |a| H(f/a) “time scaling” (XI-22) 1 |b| h(t/b) ⇐⇒ H(bf) “frequency scaling” (XI-23) h(t − t◦) ⇐⇒ H(f)e2πift◦ “time shifting” (XI-24) h(t)e−2πif◦t ⇐⇒ H(f − f◦) “freq. shifting.” (XI-25) Donald G. Luttermoser, ETSU XI–9 D. Fourier Transform of Discretely Sampled Data. 1. In the most common situations, function h(t) is sampled (i.e., measurements taken) at evenly spaced intervals in time. a) Let τ denote the time interval between consecutive sam- ples such that hn = h(nτ ) , n = . . . ,−3,−2,−1, 0, 1, 2, 3, . . . (XI-26) Often, τ is called the sampling interval. b) The reciprocal of the time interval is called the sampling rate. If τ is measured in seconds, then the sampling rate is measured in Hz (cycles per second). 2. Sampling Theorem and Aliasing. a) For any sampling interval, there is a special frequency, fc, called the Nyquist critical frequency, given by fc ≡ 1 2τ . (XI-27) i) If a sine wave of the Nyquist critical frequency is sampled at its positive peak value, then the next sample will be at the negative trough value, the sample after that at the positive peak again, and so on. ii) Expressed otherwise: Critical sampling of a sine wave is two sample points per cycle. iii) One frequently chooses to measure time in units of the sampling interval τ . In this case the Nyquist critical frequency is just the constant 1/2. XI–10 PHYS-4007/5007: Computational Physics b) The Nyquist critical frequency is important for two dis- tinct reasons, the first good, the second bad. i) The sampling theorem: If a continuous function h(t), sampled at interval τ , happens to be band- width limited to frequencies smaller than fc, then the function h(t) is completely determined by its sample hn. Explicitly h(t) = τ +∞∑ n=−∞ hn sin[2πfc(t − nτ )] π(t − nτ ) . (XI-28) ii) Sampling a continuous function that is not band- width limited to less than the Nyquist critical fre- quency will miss information outside then range of −fc < f < fc =⇒ this is called aliasing. Any fre- quency component outside of the frequency range (−fc, fc) is aliased (i.e., falsely translated) into that range by the very act of discrete sampling. The ef- fects of this are shown in Figure (XI-1). 3. Discrete Fourier Transform. a) Suppose we have N consecutive sampled values hk ≡ h(tk) , tk ≡ kτ , k = 0, 1, 2, . . . , N − 1 . (XI-29) For description purposes, let’s assume that h(t) is an even function. b) Now determine estimates of the frequency from −fc to +fc at the distinct points defined by fn ≡ n Nτ , n = −N 2 , . . . , N 2 . (XI-30) The extreme values of n correspond to the lower and upper limits of the Nyquist critical frequency range. Donald G. Luttermoser, ETSU XI–13 e) The discrete inverse Fourier transform then takes the form hk = 1 N N−1∑ n=0 Hne −2πikn/N . (XI-34) f) Finally, in analogy to Eq. (XI-16), we can write the dis- crete form of Parseval’s theorem as N−1∑ k=0 |hk|2 = 1 N N−1∑ n=0 |Hn|2. (XI-35) E. Fast Fourier Transform (FFT). 1. We can ask the question, how much time is required to compute a function of the form Hn = N−1∑ k=0 W nkhk , (XI-36) where the vector of hk’s is multiplied by a matrix whose (n, k) th is the constant W to the power n × k with W given by W ≡ e2πi/N , (XI-37) hence a Fourier series expression. a) This matrix multiplication requires N2 complex multipli- cations plus a smaller number of operations to generate the required powers of W . b) As such, discrete Fourier transforms appear then to be a O(N2) process. 2. We can speed the calculations up to order O(N log2 N) opera- tions with an algorithm called the Fast Fourier Transform or (FFT) for short. a) The difference between N log2 N and N 2 calculations is immense. XI–14 PHYS-4007/5007: Computational Physics b) With N = 106, this corresponds to 0.03 seconds and 20 minutes on a 1 GHz processor, respectively. 3. FFTs were first developed for computational coding in the mid- 1960s by J.W. Cooley and J.W. Tukey. The earliest “discoveries” of the FFT was made by Danielson and Lanczos in 1942. The Danielson-Lanczos Lemma is as follows: a) A discrete Fourier transform of length N can be rewrit- ten as the sum of two discrete Fourier transforms each of length N/2. b) One is formed from the even-numbered points of the orig- inal N , the other from the odd-numbered points. The mathematical proof is Fk = N−1∑ j=0 e2πijk/Nfj = N/2−1∑ j=0 e2πik(2j)/Nf2j + N/2−1∑ j=0 e2πik(2j+1)/Nf2j+1 = N/2−1∑ j=0 e2πikj)/(N/2)f2j + W k N/2−1∑ j=0 e2πikj/(N/2)f2j+1 = F ek + W kF ok . (XI-38) c) Note that k in the equations above varies from 0 to n, not just to N/2. Never the less, the transforms F ek (the ‘even’ sum) and F ok (the ‘odd’ sum) are periodic in k with length N/2. As such, each is repeated through 2 cycles to obtain Fk d) For this to be the most effective, N should be an integer multiple of 2. If it is not, one should pad the vectors with zeros until the next power of 2 is reached. Donald G. Luttermoser, ETSU XI–15 4. Virtually all mathematics software packages have FFTs built into the command structure or available in a math library. FFTs are most often used in the convolution and deconvolution of spectral data.