18.085 Summer 2020 Week 5 Lecture Notes

18.085/18.0851 Computational Science and Engineering I Summer 2020 Week 5 (July 6th-July 10th) Lecturer: Richard Zhang 5.1 Second Order Differential Equations: RLC Circuit, Damped Oscillation, and More The previous lecture taps into graph models and Kirkhoff laws. Those are frameworks and natural laws that eventually would boil down to a set of differential equations. More often than not, these are second-order differential equations. This lecture explores these equations in greater depth. 5.1.1 Damped Harmonic Oscillation We again start with one mass m hanging from one spring with Hooke’s constant c. The top is fixed, while the bottom is free. We let u(t) be the displacement of the mass away from equilibrium. This time, however, we assume that the system experiences some friction, proportional to the speed of the mass. In other words, f riction = ν du dt (5.1) In that case, we can write down Newton’s second law, with the newly added frictional force m d2 u du +ν + cu(t) = 0 dt dt2 u(0) = 0 u0 (0) = 1 (5.2) (5.3) (5.4) We first take the Ansatz approach as before, assuming that u(t) takes the form of : u(t) = Aebt (5.5) Before plugging the Ansatz into the differential equation, we do some prep work of computing its first and second deriatives du = bAebt dt d2 u = b2 Aebt dt2 (5.6) (5.7) Then we can plug the Ansatz into the equation and write that mb2 Aebt + bνAebt + cAebt = 0 (5.8) mb2 + bν + c = 0 (5.9) where we can cancel Aebt and obtain 5-1 5-2 Lecture 5: July 6th-July 10th This is a quadratic equation that can be solved using the formula. Hence we can write √ −ν ± ν 2 − 4mc b= 2m (5.10) For various combinations of ν, m, and c, we would have different values of b, which leads to different behaviors of the damped harmonic oscillators. 5.1.1.1 Overdamping: ν 2 − 4mc > 0 In this case, we have two real values of b. Let’s call them b1 and b2 . In other words √ −ν + ν 2 − 4mc b1 = √2m −ν − ν 2 − 4mc b2 = 2m (5.11) (5.12) Then we know that the solution to differential equations would be the sum of two cases, ie. u(t) = A1 eb1 t + A2 eb2 t (5.13) This corresponds to the case when ν is large (ie. lots of friction in the system), m is small(ie. light mass), and/or c is small (ie. weak spring). When friction dominates, there is no oscillation, but only exponential decays. This is on top of the fact that both b1 and b2 are negative. 5.1.1.2 Underdamping: ν 2 − 4mc < 0 In this case, the square root gives a complex number as we are attempting to take the square root of a √ 2 negative number. Let ω = 4mc−ν . We then write down the solution as 2m ν u(t) = e− 2m t (A1 eiωt + A2 e−iωt ) (5.14) From the structure of the solution, we can see the following two things • It has an oscillatory component, with frequency ω; • The amplitude has an exponential decay of rate 5.1.2 ν 2m RLC Circuit A typical circuit involves a resistor, an inductor, and a capacitor. If we let I be the current of the circuit, C be the capacitance of the capacitor, R be the resistance of the resistor, and L be the inductance of the inductor, we can write down the following relations between the voltage across each device and their inherent properties • Voltage across a resistor: V = IR • Voltage across a capacitor: V = 1 C R • Voltage across an inductor: V = L dI dt Idt Lecture 5: July 6th-July 10th 5-3 Kirkhoff’s voltage law says that ”in any closed loop network, the total voltage around the loop is equal to the sum of all the voltage drops within the same loop.” Hence, in a circuit with a time-dependent voltage, V (t), we can write that Z dI 1 L + RI + Idt = V (t) (5.15) dt C R If we let u(t) = I(t), we then have L du 1 d2 u +R + u = V (t) dt C dt2 (5.16) If V (t) = 0, we get the same form of equation as in the harmonic oscillator, with • L is the inertial term in the system, like the mass m in the mechanical harmonic oscillator • R is the frictional term in the system, like the frictional coefficient ν in the mechanical harmonic oscillator • 1 C is the energy storing term in the system, like the spring constant c in the mechanical harmonic oscillator 5.2 The Laplace Transform The method of the Ansatz, described above, essentially borrows into the ”frequency” space before reassembling the solution in the original time space. The method can be formalized into a method called the Laplace transform, which transforms the entire problem in the frequency domain. Then in the last step, we transform it back to the time domain. Here is the definition of the Laplace transform Definition 5.1 For a given function u(t), the Laplace transform, U (s), is defined as Z ∞ U (s) = u(t)e−st dt (5.17) 0 As an example, the Laplace transform of the function, u(t) = 1, would be Z ∞ U (s) = e−st dt (5.18) 0 =− = 5.2.1 1 s e−st ∞ | s 0 (5.19) (5.20) Integration by Parts Another mathematical tool that we will need to understand and use the Laplace transform to solve differential equations would be integration by parts. For those of you who have learned integration by parts in calculus, this will be a review class. 5-4 Lecture 5: July 6th-July 10th Everything starts with the product rule. Supposed I have my functions, f (x) and g(x). Then we can write the derivatives of f ∗ g as d df dg (f g) = g+ f dx dx dx In the short-hand notation, we can write that (f g)0 = f 0 g + f g 0 (5.22) Now suppose we perform integration over [a, b] on both sides, where −∞ ≤ a < b ≤ ∞. Then Z a Z a a 0 f g|b = f gdx + f g 0 dx b (5.21) (5.23) b An algebraic step would give us a Z 0 f gdx = b f g|ab a Z f g 0 dx − (5.24) b Here is an application of integration by parts. Suppose I want to integrate ex sin x. We do one step of integration by parts, assuming that f 0 = ex and g = sin(x). Then according to the formula, Z a Z ex sin x = ex sin x − ex (− cos x)dx (5.25) b Z = ex sin x + ex cos xdx (5.26) Now we perform another round of integration by parts, assuming f 0 = ex and g = cos x. Then Z a Z x x x e sin x = e sin x + e cos(x) − ex sin xdx b Z a 2 ex sin x = ex sin x + ex cos x Zb a 1 ex sin x = (ex sin x + ex cos x) 2 b (5.27) (5.28) (5.29) Now we are ready to study the Laplace transform 5.2.2 Application of the Laplace Transform to Solving Differential Equations Let’s see how we use the Laplace transform to solve the differential equation m 5.2.2.1 d2 u du +ν + cu = 0 2 dt dt u(0) = 0 u0 (0) = 1 (5.30) (5.31) (5.32) Step 1: Apply the Laplace transform on both sides R∞ We multiply both sides by e−st and apply the integral 0 dt Z ∞ u Z ∞ Z ∞ du −st d −st m e dt + ν e dt + c ue−st dt = 0 2 dt dt 0 0 0 (5.33) (5.34) Lecture 5: July 6th-July 10th Since by definition, R∞ 0 5-5 u(t)e−st dt, we perform the integration by parts on the friction term and obtain Z ∞ Z ∞ du −st −st ∞ ue−st dt (5.35) e dt = ue |0 + s dt 0 0 = −u(0) + sU (s) (5.36) = sU (s) (5.37) Then then we perform the integration by parts twice on the inertial term and obtain that Z ∞ 2 Z ∞ d u −st du −st ∞ du −st e dt = e | + s e dt 0 dt dt dt2 0 0 Z ∞ du(0) =− + su(0) + s2 ue−st dt dt 0 = −1 + s2 U (s) 5.2.2.2 (5.38) (5.39) (5.40) Step 2: Write down the Laplace transformed equation and massage the algebra if needed Putting everything above together, we obtain that −m + ms2 U (s) + νsU (s) + cU (s) = 0 (5.41) 2 (ms + νs + c)U (s) − m = 0 (5.42) m U (s) = 2 ms + νs + c 1 U (s) = 2 s + (ν/m)s + c/m (5.43) (5.44) For the purpose that will be obvious later, we shall do a partial fraction decomposition on U (s). Let a, b ∈ C be the complex roots of s2 + (ν/m)s + c/m. Then U (s) can be expressed as U (s) = 1 (s − a)(s − b) (5.45) The partial fraction decomposition assumes that U (s) can be expressed as U (s) = A B + s−a s−b (5.46) where A and B are constants, such that A(s − b) + B(s − a) = 1 (5.47) A+B =0 Ab + Ba = −1 (5.48) (5.49) whence, from which we can see that B = −A and A = U (s) = 1 a−b . 1 a−b Hence we write  1 1 − s−a s−b  (5.50) 5-6 5.2.2.3 Lecture 5: July 6th-July 10th Step 3: ”Perform” the inverse Laplace transform Now we just have to transform U (s) back to u(t). The actual operation of the inverse Laplace transform is difficult, and we will not get into the details for this course. There are mathematicians who have worked out the Laplace transforms and their inverses and curated the results into a table. We just have to use it. 1 1 Here we need the inverse Laplace transform of s−a and s−b . According to the table, the inverse Laplace a b transform are e and e respectively. Hence we can write that u(t) = 1 (eat − ebt ) a−b (5.51) which is what we would get upon solving the initial value problem. This will be left as a homework problem. 5.2.3 A note about nondimensionalizing equations The equation above is dimensional, which is slight awkward to analyze because of the units floating around. A good engineer would know the importance of nondimensionalizing equations immediately. here is how you would nondimensionalize the damped harmonic oscillator m d2 u du +ν + cu = 0 dt dt2 (5.52) (5.53) First we divide both sides by m ν du c d2 u + + u=0 2 m dt m dt (5.54) (5.55) Then we fix some characteristic length l and we let q = u/l. Then the equation becomes ν dq c d2 q + + q=0 m dt m dt2 (5.56) (5.57) We let γ = ν 2m and ω0 = pc m. Finally the equation becomes d2 q dq + 2γ + ω02 q = 0 dt dt2 (5.58) (5.59) Note that γ and ω0 have the unit of frequency. There will be a homework problem on repeating the same analysis of damped harmonic oscillator using this nondimensional equation. 5.2.4 Another Example: Driven Oscillator So far we’ve been dealing with unforced oscillator. Let’s examine the driven oscillator. Suppose we have the following damped harmonic oscillator u00 + 2.05u0 + u = 1 u(0) = 0 u0 (0) = 0 We can solve this using the Laplace transform (5.60) (5.61) (5.62) Lecture 5: July 6th-July 10th • U (s) = R∞ 0 5-7 e−st u(t)dt • sU (s) − u(0) = sU (s) = R∞ • s2 U (s) − su(0) − u0 (0) = 0 u0 (t)e−st ds R∞ 0 u00 (t)e−st ds Since u(0) = u0 (0) = 0, we have s2 U (s) + 2.05sU (s) + U (s) = 1/s   1 1 U (s) = s s2 + 2.05s + 1 In engineering, 1 s2 +2.05s+1 (5.63) (5.64) is called the transfer function. It connects the input with the output. To transform it back, we need to perform a partial fraction decomposition on 1 s(s + 0.8)(s + 1.25) 1 −25/9 16/9 = + + s s + 0.8 s + 1.25 U (s) = (5.65) (5.66) Hence, u(t) = 1 − (25/9)e−0.8t + (16/9)e−1.25t (5.67) A few comments • Without the forcing term 1, this is a purely overdamped system whose amplitude decays to zero as t goes to infinity • With the forcing term, however, the system starts from zero and is growing to one, where it stays • 1 is therefore the steady state solution, while the decaying exponential is the transient state solution. This is going to be true in general: the solution to a driven oscillator is the sum of transient state plus steady state solutions.