<jupyter_start><jupyter_text>PHY321: Harmonic Oscillations, Damping, Resonances and time-dependent Forces  **[<NAME>](http://mhjgit.github.io/info/doc/web/)**, Department of Physics and Astronomy and Facility for Rare Ion Beams (FRIB), Michigan State University, USA and Department of Physics, University of Oslo, NorwayDate: **Mar 1, 2021**Copyright 1999-2021, [<NAME>](http://mhjgit.github.io/info/doc/web/). Released under CC Attribution-NonCommercial 4.0 license Aims and Overarching Motivation MondayDamped oscillations. Analytical and numerical solutions**Reading suggestion**: Taylor sections 5.4-5.5.  WednesdayNo lecture, study day FridayDriven oscillations and resonances with examples.**Reading suggestion**: Taylor sections 5.5-5.6.  Damped OscillatorsWe consider only the case where the damping force is proportional tothe velocity. This is counter to dragging friction, where the force isproportional in strength to the normal force and independent ofvelocity, and is also inconsistent with wind resistance, where themagnitude of the drag force is proportional the square of thevelocity. Rolling resistance does seem to be mainly proportional tothe velocity. However, the main motivation for considering dampingforces proportional to the velocity is that the math is morefriendly. This is because the differential equation is linear,i.e. each term is of order $x$, $\dot{x}$, $\ddot{x}\cdots$, or eventerms with no mention of $x$, and there are no terms such as $x^2$ or$x\ddot{x}$. The equations of motion for a spring with damping force$-b\dot{x}$ are $$\begin{equation}m\ddot{x}+b\dot{x}+kx=0.\label{_auto1} \tag{1}\end{equation}$$  Harmonic Oscillator,  DampingJust to make the solution a bit less messy, we rewrite this equation as $$\begin{equation}\label{eq:dampeddiffyq} \tag{2}\ddot{x}+2\beta\dot{x}+\omega_0^2x=0,~~~~\beta\equiv b/2m,~\omega_0\equiv\sqrt{k/m}.\end{equation}$$ Both $\beta$ and $\omega$ have dimensions of inverse time. To find solutions (see appendix C in the text) you must make an educated guess at the form of the solution. To do this, first realize that the solution will need an arbitrary normalization $A$ because the equation is linear. Secondly, realize that if the form is $$\begin{equation}x=Ae^{rt}\label{_auto2} \tag{3}\end{equation}$$ that each derivative simply brings out an extra power of $r$. Thismeans that the $Ae^{rt}$ factors out and one can simply solve for anequation for $r$. Plugging this form into Eq. ([2](eq:dampeddiffyq)), $$\begin{equation}r^2+2\beta r+\omega_0^2=0.\label{_auto3} \tag{4}\end{equation}$$  Harmonic Oscillator,  Solutions of Damped MotionBecause this is a quadratic equation there will be two solutions, $$\begin{equation}r=-\beta\pm\sqrt{\beta^2-\omega_0^2}.\label{_auto4} \tag{5}\end{equation}$$ We refer to the two solutions as $r_1$ and $r_2$ corresponding to the$+$ and $-$ roots. As expected, there should be two arbitraryconstants involved in the solution, $$\begin{equation}x=A_1e^{r_1t}+A_2e^{r_2t},\label{_auto5} \tag{6}\end{equation}$$ where the coefficients $A_1$ and $A_2$ are determined by initialconditions.The roots listed above, $\sqrt{\omega_0^2-\beta_0^2}$, will beimaginary if the damping is small and $\beta<\omega_0$. In that case,$r$ is complex and the factor $\exp{(rt)}$ will have some oscillatorybehavior. If the roots are real, there will only be exponentiallydecaying solutions. There are three cases: Underdamped: $\beta<\omega_0$ $$\begin{eqnarray}x&=&A_1e^{-\beta t}e^{i\omega't}+A_2e^{-\beta t}e^{-i\omega't},~~\omega'\equiv\sqrt{\omega_0^2-\beta^2}\\\nonumber&=&(A_1+A_2)e^{-\beta t}\cos\omega't+i(A_1-A_2)e^{-\beta t}\sin\omega't.\end{eqnarray}$$ Here we have made use of the identity$e^{i\omega't}=\cos\omega't+i\sin\omega't$. Because the constants arearbitrary, and because the real and imaginary parts are both solutionsindividually, we can simply consider the real part of the solutionalone: $$\begin{eqnarray}\label{eq:homogsolution} \tag{7}x&=&B_1e^{-\beta t}\cos\omega't+B_2e^{-\beta t}\sin\omega't,\\\nonumber \omega'&\equiv&\sqrt{\omega_0^2-\beta^2}.\end{eqnarray}$$  Critical dampling: $\beta=\omega_0$In this case the two terms involving $r_1$ and $r_2$ are identicalbecause $\omega'=0$. Because we need to arbitrary constants, thereneeds to be another solution. This is found by simply guessing, or bytaking the limit of $\omega'\rightarrow 0$ from the underdampedsolution. The solution is then $$\begin{equation}\label{eq:criticallydamped} \tag{8}x=Ae^{-\beta t}+Bte^{-\beta t}.\end{equation}$$ The critically damped solution is interesting because the solutionapproaches zero quickly, but does not oscillate. For a problem withzero initial velocity, the solution never crosses zero. This is a goodchoice for designing shock absorbers or swinging doors. Overdamped: $\beta>\omega_0$ $$\begin{eqnarray}x&=&A_1\exp{-(\beta+\sqrt{\beta^2-\omega_0^2})t}+A_2\exp{-(\beta-\sqrt{\beta^2-\omega_0^2})t}\end{eqnarray}$$ This solution will also never pass the origin more than once, and thenonly if the initial velocity is strong and initially toward zero.Given $b$, $m$ and $\omega_0$, find $x(t)$ for a particle whoseinitial position is $x=0$ and has initial velocity $v_0$ (assuming anunderdamped solution).The solution is of the form, $$\begin{eqnarray*}x&=&e^{-\beta t}\left[A_1\cos(\omega' t)+A_2\sin\omega't\right],\\\dot{x}&=&-\beta x+\omega'e^{-\beta t}\left[-A_1\sin\omega't+A_2\cos\omega't\right].\\\omega'&\equiv&\sqrt{\omega_0^2-\beta^2},~~~\beta\equiv b/2m.\end{eqnarray*}$$ From the initial conditions, $A_1=0$ because $x(0)=0$ and $\omega'A_2=v_0$. So $$x=\frac{v_0}{\omega'}e^{-\beta t}\sin\omega't.$$  Harmonic Oscillator,  SolutionsConsider a single solution with no arbitrary constants, which we willcall a **particular solution**, $x_p(t)$. It should be emphasizedthat this is **A** particular solution, because there exists aninfinite number of such solutions because the general solution shouldhave two arbitrary constants. Now consider solutions to the sameequation without the driving term, which include two arbitraryconstants. These are called either **homogenous solutions** or **complementary solutions**, and were given in the previous section,e.g. Eq. ([7](eq:homogsolution)) for the underdamped case. Thehomogenous solution already incorporates the two arbitrary constants,so any sum of a homogenous solution and a particular solution willrepresent the **general solution** of the equation. The generalsolution incorporates the two arbitrary constants $A$ and $B$ toaccommodate the two initial conditions. One could have picked adifferent particular solution, i.e. the original particular solutionplus any homogenous solution with the arbitrary constants $A_p$ and$B_p$ chosen at will. When one adds in the homogenous solution, whichhas adjustable constants with arbitrary constants $A'$ and $B'$, tothe new particular solution, one can get the same general solution bysimply adjusting the new constants such that $A'+A_p=A$ and$B'+B_p=B$. Thus, the choice of $A_p$ and $B_p$ are irrelevant, andwhen choosing the particular solution it is best to make the simplestchoice possible. Harmonic Oscillator,  Particular SolutionTo find a particular solution, one first guesses at the form, $$\begin{equation}\label{eq:partform} \tag{9}x_p(t)=D\cos(\omega t-\delta),\end{equation}$$ and rewrite the differential equation as $$\begin{equation}D\left\{-\omega^2\cos(\omega t-\delta)-2\beta\omega\sin(\omega t-\delta)+\omega_0^2\cos(\omega t-\delta)\right\}=\frac{F_0}{m}\cos(\omega t).\label{_auto6} \tag{10}\end{equation}$$ One can now use angle addition formulas to get $$\begin{eqnarray}D\left\{(-\omega^2\cos\delta+2\beta\omega\sin\delta+\omega_0^2\cos\delta)\cos(\omega t)\right.&&\\\nonumber\left.+(-\omega^2\sin\delta-2\beta\omega\cos\delta+\omega_0^2\sin\delta)\sin(\omega t)\right\}&=&\frac{F_0}{m}\cos(\omega t).\end{eqnarray}$$ Both the $\cos$ and $\sin$ terms need to equate if the expression is to hold at all times. Thus, this becomes two equations $$\begin{eqnarray}D\left\{-\omega^2\cos\delta+2\beta\omega\sin\delta+\omega_0^2\cos\delta\right\}&=&\frac{F_0}{m}\\\nonumber-\omega^2\sin\delta-2\beta\omega\cos\delta+\omega_0^2\sin\delta&=&0.\end{eqnarray}$$ After dividing by $\cos\delta$, the lower expression leads to $$\begin{equation}\tan\delta=\frac{2\beta\omega}{\omega_0^2-\omega^2}.\label{_auto7} \tag{11}\end{equation}$$  Solving with Driven OscillationsUsing the identities $\tan^2+1=\csc^2$ and $\sin^2+\cos^2=1$, one can also express $\sin\delta$ and $\cos\delta$, $$\begin{eqnarray}\sin\delta&=&\frac{2\beta\omega}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}},\\\nonumber\cos\delta&=&\frac{(\omega_0^2-\omega^2)}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}}\end{eqnarray}$$ Inserting the expressions for $\cos\delta$ and $\sin\delta$ into the expression for $D$, $$\begin{equation}\label{eq:Ddrive} \tag{12}D=\frac{F_0/m}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}}.\end{equation}$$ For a given initial condition, e.g. initial displacement and velocity,one must add the homogenous solution then solve for the two arbitraryconstants. However, because the homogenous solutions decay with timeas $e^{-\beta t}$, the particular solution is all that remains atlarge times, and is therefore the steady state solution. Because thearbitrary constants are all in the homogenous solution, all memory ofthe initial conditions are lost at large times, $t>>1/\beta$.The amplitude of the motion, $D$, is linearly proportional to thedriving force ($F_0/m$), but also depends on the driving frequency$\omega$. For small $\beta$ the maximum will occur at$\omega=\omega_0$. This is referred to as a resonance. In the limit$\beta\rightarrow 0$ the amplitude at resonance approaches infinity. Alternative Derivation for Driven OscillatorsHere, we derive the same expressions as in Equations ([9](eq:partform)) and ([12](eq:Ddrive)) but express the driving forces as $$\begin{eqnarray}F(t)&=&F_0e^{i\omega t},\end{eqnarray}$$ rather than as $F_0\cos\omega t$. The real part of $F$ is the same as before. For the differential equation, $$\begin{eqnarray}\label{eq:compdrive} \tag{13}\ddot{x}+2\beta\dot{x}+\omega_0^2x&=&\frac{F_0}{m}e^{i\omega t},\end{eqnarray}$$ one can treat $x(t)$ as an imaginary function. Because the operations$d^2/dt^2$ and $d/dt$ are real and thus do not mix the real andimaginary parts of $x(t)$, Eq. ([13](eq:compdrive)) is effectively 2equations. Because $e^{\omega t}=\cos\omega t+i\sin\omega t$, the realpart of the solution for $x(t)$ gives the solution for a driving force$F_0\cos\omega t$, and the imaginary part of $x$ corresponds to thecase where the driving force is $F_0\sin\omega t$. It is rather easyto solve for the complex $x$ in this case, and by taking the real partof the solution, one finds the answer for the $\cos\omega t$ drivingforce.We assume a simple form for the particular solution $$\begin{equation}x_p=De^{i\omega t},\label{_auto8} \tag{14}\end{equation}$$ where $D$ is a complex constant.From Eq. ([13](eq:compdrive)) one inserts the form for $x_p$ above to get $$\begin{eqnarray}D\left\{-\omega^2+2i\beta\omega+\omega_0^2\right\}e^{i\omega t}=(F_0/m)e^{i\omega t},\\\nonumberD=\frac{F_0/m}{(\omega_0^2-\omega^2)+2i\beta\omega}.\end{eqnarray}$$ The norm and phase for $D=|D|e^{-i\delta}$ can be read by inspection, $$\begin{equation}|D|=\frac{F_0/m}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}},~~~~\tan\delta=\frac{2\beta\omega}{\omega_0^2-\omega^2}.\label{_auto9} \tag{15}\end{equation}$$ This is the same expression for $\delta$ as before. One then finds $x_p(t)$, $$\begin{eqnarray}\label{eq:fastdriven1} \tag{16}x_p(t)&=&\Re\frac{(F_0/m)e^{i\omega t-i\delta}}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}\\\nonumber&=&\frac{(F_0/m)\cos(\omega t-\delta)}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}.\end{eqnarray}$$ This is the same answer as before.If one wished to solve for the case where $F(t)= F_0\sin\omega t$, the imaginary part of the solution would work $$\begin{eqnarray}\label{eq:fastdriven2} \tag{17}x_p(t)&=&\Im\frac{(F_0/m)e^{i\omega t-i\delta}}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}\\\nonumber&=&\frac{(F_0/m)\sin(\omega t-\delta)}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}.\end{eqnarray}$$  Damped and Driven OscillatorConsider the damped and driven harmonic oscillator worked out above. Given $F_0, m,\beta$ and $\omega_0$, solve for the complete solution $x(t)$ for the case where $F=F_0\sin\omega t$ with initial conditions $x(t=0)=0$ and $v(t=0)=0$. Assume the underdamped case.The general solution including the arbitrary constants includes both the homogenous and particular solutions, $$\begin{eqnarray*}x(t)&=&\frac{F_0}{m}\frac{\sin(\omega t-\delta)}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}+A\cos\omega't e^{-\beta t}+B\sin\omega't e^{-\beta t}.\end{eqnarray*}$$ The quantities $\delta$ and $\omega'$ are given earlier in thesection, $\omega'=\sqrt{\omega_0^2-\beta^2},\delta=\tan^{-1}(2\beta\omega/(\omega_0^2-\omega^2)$. Here, solvingthe problem means finding the arbitrary constants $A$ and$B$. Satisfying the initial conditions for the initial position andvelocity: $$\begin{eqnarray*}x(t=0)=0&=&-\eta\sin\delta+A,\\v(t=0)=0&=&\omega\eta\cos\delta-\beta A+\omega'B,\\\eta&\equiv&\frac{F_0}{m}\frac{1}{\sqrt{(\omega_0^2-\omega^2)^2+4\beta^2\omega^2}}.\end{eqnarray*}$$ The problem is now reduced to 2 equations and 2 unknowns, $A$ and $B$. The solution is $$\begin{eqnarray}A&=& \eta\sin\delta ,~~~B=\frac{-\omega\eta\cos\delta+\beta\eta\sin\delta}{\omega'}.\end{eqnarray}$$  Resonance Widths; the $Q$ factorFrom the previous two sections, the particular solution for a driving force, $F=F_0\cos\omega t$, is $$\begin{eqnarray}x_p(t)&=&\frac{F_0/m}{\sqrt{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}}\cos(\omega_t-\delta),\\\nonumber\delta&=&\tan^{-1}\left(\frac{2\beta\omega}{\omega_0^2-\omega^2}\right).\end{eqnarray}$$ If one fixes the driving frequency $\omega$ and adjusts thefundamental frequency $\omega_0=\sqrt{k/m}$, the maximum amplitudeoccurs when $\omega_0=\omega$ because that is when the term from thedenominator $(\omega_0^2-\omega^2)^2+4\omega^2\beta^2$ is at aminimum. This is akin to dialing into a radio station. However, if onefixes $\omega_0$ and adjusts the driving frequency one minimize withrespect to $\omega$, e.g. set $$\begin{equation}\frac{d}{d\omega}\left[(\omega_0^2-\omega^2)^2+4\omega^2\beta^2\right]=0,\label{_auto10} \tag{18}\end{equation}$$ and one finds that the maximum amplitude occurs when$\omega=\sqrt{\omega_0^2-2\beta^2}$. If $\beta$ is small relative to$\omega_0$, one can simply state that the maximum amplitude is $$\begin{equation}x_{\rm max}\approx\frac{F_0}{2m\beta \omega_0}.\label{_auto11} \tag{19}\end{equation}$$ $$\begin{eqnarray}\frac{4\omega^2\beta^2}{(\omega_0^2-\omega^2)^2+4\omega^2\beta^2}=\frac{1}{2}.\end{eqnarray}$$ For small damping this occurs when $\omega=\omega_0\pm \beta$, so the $FWHM\approx 2\beta$. For the purposes of tuning to a specific frequency, one wants the width to be as small as possible. The ratio of $\omega_0$ to $FWHM$ is known as the _quality_factor, or $Q$ factor, $$\begin{equation}Q\equiv \frac{\omega_0}{2\beta}.\label{_auto12} \tag{20}\end{equation}$$  Numerical Studies of Driven OscillationsSolving the problem of driven oscillations numerically gives us muchmore flexibility to study different types of driving forces. We canreuse our earlier code by simply adding a driving force. If we stay inthe $x$-direction only this can be easily done by adding a term$F_{\mathrm{ext}}(x,t)$. Note that we have kept it rather generalhere, allowing for both a spatial and a temporal dependence.Before we dive into the code, we need to briefly remind ourselvesabout the equations we started with for the case with damping, namely $$m\frac{d^2x}{dt^2} + b\frac{dx}{dt}+kx(t) =0,$$ with no external force applied to the system.Let us now for simplicty assume that our external force is given by $$F_{\mathrm{ext}}(t) = F_0\cos{(\omega t)},$$ where $F_0$ is a constant (what is its dimension?) and $\omega$ is the frequency of the applied external driving force.**Small question:** would you expect energy to be conserved now?Introducing the external force into our lovely differential equationand dividing by $m$ and introducing $\omega_0^2=\sqrt{k/m}$ we have $$\frac{d^2x}{dt^2} + \frac{b}{m}\frac{dx}{dt}+\omega_0^2x(t) =\frac{F_0}{m}\cos{(\omega t)},$$ Thereafter we introduce a dimensionless time $\tau = t\omega_0$and a dimensionless frequency $\tilde{\omega}=\omega/\omega_0$. We have then $$\frac{d^2x}{d\tau^2} + \frac{b}{m\omega_0}\frac{dx}{d\tau}+x(\tau) =\frac{F_0}{m\omega_0^2}\cos{(\tilde{\omega}\tau)},$$ Introducing a new amplitude $\tilde{F} =F_0/(m\omega_0^2)$ (check dimensionality again) we have $$\frac{d^2x}{d\tau^2} + \frac{b}{m\omega_0}\frac{dx}{d\tau}+x(\tau) =\tilde{F}\cos{(\tilde{\omega}\tau)}.$$ Our final step, as we did in the case of various types of damping, isto define $\gamma = b/(2m\omega_0)$ and rewrite our equations as $$\frac{d^2x}{d\tau^2} + 2\gamma\frac{dx}{d\tau}+x(\tau) =\tilde{F}\cos{(\tilde{\omega}\tau)}.$$ This is the equation we will code below using the Euler-Cromer method.<jupyter_code>DeltaT = 0.001
#set up arrays 
tfinal = 20 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, v, and x
t = np.zeros(n)
v = np.zeros(n)
x = np.zeros(n)
# Initial conditions as one-dimensional arrays of time
x0 =  1.0 
v0 = 0.0
x[0] = x0
v[0] = v0
gamma = 0.2
Omegatilde = 0.5
Ftilde = 1.0
# Start integrating using Euler-Cromer's method
for i in range(n-1):
    # Set up the acceleration
    # Here you could have defined your own function for this
    a =  -2*gamma*v[i]-x[i]+Ftilde*cos(t[i]*Omegatilde)
    # update velocity, time and position
    v[i+1] = v[i] + DeltaT*a
    x[i+1] = x[i] + DeltaT*v[i+1]
    t[i+1] = t[i] + DeltaT
# Plot position as function of time    
fig, ax = plt.subplots()
ax.set_ylabel('x[m]')
ax.set_xlabel('t[s]')
ax.plot(t, x)
fig.tight_layout()
save_fig("ForcedBlockEulerCromer")
plt.show()<jupyter_output><empty_output><jupyter_text>In the above example we have focused on the Euler-Cromer method. Thismethod has a local truncation error which is proportional to $\Delta t^2$and thereby a global error which is proportional to $\Delta t$.We can improve this by using the Runge-Kutta family ofmethods. The widely popular Runge-Kutta to fourth order or just **RK4**has indeed a much better truncation error. The RK4 method has a globalerror which is proportional to $\Delta t$.Let us revisit this method and see how we can implement it for the above example. Differential Equations, Runge-Kutta methodsRunge-Kutta (RK) methods are based on Taylor expansion formulae, but yieldin general better algorithms for solutions of an ordinary differential equation.The basic philosophy is that it provides an intermediate step in the computation of $y_{i+1}$.To see this, consider first the following definitions $$\begin{equation}\frac{dy}{dt}=f(t,y),  \label{_auto13} \tag{21}\end{equation}$$ and $$\begin{equation}y(t)=\int f(t,y) dt,  \label{_auto14} \tag{22}\end{equation}$$ and $$\begin{equation}y_{i+1}=y_i+ \int_{t_i}^{t_{i+1}} f(t,y) dt.\label{_auto15} \tag{23}\end{equation}$$ To demonstrate the philosophy behind RK methods, let us considerthe second-order RK method, RK2.The first approximation consists in Taylor expanding $f(t,y)$around the center of the integration interval $t_i$ to $t_{i+1}$,that is, at $t_i+h/2$, $h$ being the step.Using the midpoint formula for an integral, defining $y(t_i+h/2) = y_{i+1/2}$ and   $t_i+h/2 = t_{i+1/2}$, we obtain $$\begin{equation}\int_{t_i}^{t_{i+1}} f(t,y) dt \approx hf(t_{i+1/2},y_{i+1/2}) +O(h^3).\label{_auto16} \tag{24}\end{equation}$$ This means in turn that we have $$\begin{equation}y_{i+1}=y_i + hf(t_{i+1/2},y_{i+1/2}) +O(h^3).\label{_auto17} \tag{25}\end{equation}$$ However, we do not know the value of   $y_{i+1/2}$. Here comes thus the next approximation, namely, we use Euler'smethod to approximate $y_{i+1/2}$. We have then $$\begin{equation}y_{(i+1/2)}=y_i + \frac{h}{2}\frac{dy}{dt}=y(t_i) + \frac{h}{2}f(t_i,y_i).\label{_auto18} \tag{26}\end{equation}$$ This means that we can define the following algorithm for the second-order Runge-Kutta method, RK2. 46 <<<!!MATH_BLOCK $$\begin{equation}k_2=hf(t_{i+1/2},y_i+k_1/2),\label{_auto20} \tag{28}\end{equation}$$ with the final value $$\begin{equation} y_{i+i}\approx y_i + k_2 +O(h^3). \label{_auto21} \tag{29}\end{equation}$$ The difference between the previous one-step methods is that we now need an intermediate step in our evaluation,namely $t_i+h/2 = t_{(i+1/2)}$ where we evaluate the derivative $f$. This involves more operations, but the gain is a better stabilityin the solution.The fourth-order Runge-Kutta, RK4, has the following algorithm 49 <<<!!MATH_BLOCK $$k_3=hf(t_i+h/2,y_i+k_2/2)\hspace{0.5cm}   k_4=hf(t_i+h,y_i+k_3)$$ with the final result $$y_{i+1}=y_i +\frac{1}{6}\left( k_1 +2k_2+2k_3+k_4\right).$$ Thus, the algorithm consists in first calculating $k_1$ with $t_i$, $y_1$ and $f$ as inputs. Thereafter, we increase the stepsize by $h/2$ and calculate $k_2$, then $k_3$ and finally $k_4$. The global error goes as $O(h^4)$.However, at this stage, if we keep adding different methods in ourmain program, the code will quickly become messy and ugly. Before weproceed thus, we will now introduce functions that enbody the variousmethods for solving differential equations. This means that we canseparate out these methods in own functions and files (and later as classes and moregeneric functions) and simply call them when needed. Similarly, wecould easily encapsulate various forces or other quantities ofinterest in terms of functions. To see this, let us bring up the codewe developed above for the simple sliding block, but now only with the simple forward Euler method. We introducetwo functions, one for the simple Euler method and one for theforce.Note that here the forward Euler method does not know the specific force function to be called.It receives just an input the name. We can easily change the force by adding another function.<jupyter_code>def ForwardEuler(v,x,t,n,Force):
    for i in range(n-1):
        v[i+1] = v[i] + DeltaT*Force(v[i],x[i],t[i])
        x[i+1] = x[i] + DeltaT*v[i]
        t[i+1] = t[i] + DeltaT
def SpringForce(v,x,t):
#   note here that we have divided by mass and we return the acceleration
    return  -2*gamma*v-x+Ftilde*cos(t*Omegatilde)<jupyter_output><empty_output><jupyter_text>It is easy to add a new method like the Euler-Cromer<jupyter_code>def ForwardEulerCromer(v,x,t,n,Force):
    for i in range(n-1):
        a = Force(v[i],x[i],t[i])
        v[i+1] = v[i] + DeltaT*a
        x[i+1] = x[i] + DeltaT*v[i+1]
        t[i+1] = t[i] + DeltaT<jupyter_output><empty_output><jupyter_text>and the Velocity Verlet method (be careful with time-dependence here, it is not an ideal method for non-conservative forces))<jupyter_code>def VelocityVerlet(v,x,t,n,Force):
    for i in range(n-1):
        a = Force(v[i],x[i],t[i])
        x[i+1] = x[i] + DeltaT*v[i]+0.5*a
        anew = Force(v[i],x[i+1],t[i+1])
        v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
        t[i+1] = t[i] + DeltaT<jupyter_output><empty_output><jupyter_text>Finally, we can now add the Runge-Kutta2 method via a new function<jupyter_code>def RK2(v,x,t,n,Force):
    for i in range(n-1):
# Setting up k1
        k1x = DeltaT*v[i]
        k1v = DeltaT*Force(v[i],x[i],t[i])
# Setting up k2
        vv = v[i]+k1v*0.5
        xx = x[i]+k1x*0.5
        k2x = DeltaT*vv
        k2v = DeltaT*Force(vv,xx,t[i]+DeltaT*0.5)
# Final result
        x[i+1] = x[i]+k2x
        v[i+1] = v[i]+k2v
	t[i+1] = t[i]+DeltaT<jupyter_output><empty_output><jupyter_text>Finally, we can now add the Runge-Kutta2 method via a new function<jupyter_code>def RK4(v,x,t,n,Force):
    for i in range(n-1):
# Setting up k1
        k1x = DeltaT*v[i]
        k1v = DeltaT*Force(v[i],x[i],t[i])
# Setting up k2
        vv = v[i]+k1v*0.5
        xx = x[i]+k1x*0.5
        k2x = DeltaT*vv
        k2v = DeltaT*Force(vv,xx,t[i]+DeltaT*0.5)
# Setting up k3
        vv = v[i]+k2v*0.5
        xx = x[i]+k2x*0.5
        k3x = DeltaT*vv
        k3v = DeltaT*Force(vv,xx,t[i]+DeltaT*0.5)
# Setting up k4
        vv = v[i]+k3v
        xx = x[i]+k3x
        k4x = DeltaT*vv
        k4v = DeltaT*Force(vv,xx,t[i]+DeltaT)
# Final result
        x[i+1] = x[i]+(k1x+2*k2x+2*k3x+k4x)/6.
        v[i+1] = v[i]+(k1v+2*k2v+2*k3v+k4v)/6.
        t[i+1] = t[i] + DeltaT<jupyter_output><empty_output><jupyter_text>The Runge-Kutta family of methods are particularly useful when we have a time-dependent acceleration.If we have forces which depend only the spatial degrees of freedom (no velocity and/or time-dependence), then energy conserving methods like the Velocity Verlet or the Euler-Cromer method are preferred. As soon as we introduce an explicit time-dependence and/or add dissipitave forces like friction or air resistance, then methods like the family of Runge-Kutta methods are well suited for this. The code below uses the Runge-Kutta4 methods.<jupyter_code>DeltaT = 0.001
#set up arrays 
tfinal = 20 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, v, and x
t = np.zeros(n)
v = np.zeros(n)
x = np.zeros(n)
# Initial conditions (can change to more than one dim)
x0 =  1.0 
v0 = 0.0
x[0] = x0
v[0] = v0
gamma = 0.2
Omegatilde = 0.5
Ftilde = 1.0
# Start integrating using Euler's method
# Note that we define the force function as a SpringForce
RK4(v,x,t,n,SpringForce)

# Plot position as function of time    
fig, ax = plt.subplots()
ax.set_ylabel('x[m]')
ax.set_xlabel('t[s]')
ax.plot(t, x)
fig.tight_layout()
save_fig("ForcedBlockRK4")
plt.show()<jupyter_output><empty_output><jupyter_text>Principle of Superposition and Periodic Forces (Fourier Transforms)If one has several driving forces, $F(t)=\sum_n F_n(t)$, one can findthe particular solution to each $F_n$, $x_{pn}(t)$, and the particularsolution for the entire driving force is $$\begin{equation}x_p(t)=\sum_nx_{pn}(t).\label{_auto22} \tag{30}\end{equation}$$ This is known as the principal of superposition. It only applies whenthe homogenous equation is linear. If there were an anharmonic termsuch as $x^3$ in the homogenous equation, then when one summed varioussolutions, $x=(\sum_n x_n)^2$, one would get crossterms. Superposition is especially useful when $F(t)$ can be writtenas a sum of sinusoidal terms, because the solutions for eachsinusoidal (sine or cosine)  term is analytic, as we saw above.Driving forces are often periodic, even when they are notsinusoidal. Periodicity implies that for some time $\tau$ $$\begin{eqnarray}F(t+\tau)=F(t). \end{eqnarray}$$ One example of a non-sinusoidal periodic force is a square wave. Manycomponents in electric circuits are non-linear, e.g. diodes, whichmakes many wave forms non-sinusoidal even when the circuits are beingdriven by purely sinusoidal sources.The code here shows a typical example of such a square wave generated using the functionality included in the **scipy** Python package. We have used a period of $\tau=0.2$.<jupyter_code>%matplotlib inline

import numpy as np
import math
from scipy import signal
import matplotlib.pyplot as plt

# number of points                                                                                       
n = 500
# start and final times                                                                                  
t0 = 0.0
tn = 1.0
# Period                                                                                                 
t = np.linspace(t0, tn, n, endpoint=False)
SqrSignal = np.zeros(n)
SqrSignal = 1.0+signal.square(2*np.pi*5*t)
plt.plot(t, SqrSignal)
plt.ylim(-0.5, 2.5)
plt.show()<jupyter_output><empty_output><jupyter_text>For the sinusoidal example studied in the previous subsections theperiod is $\tau=2\pi/\omega$. However, higher harmonics can alsosatisfy the periodicity requirement. In general, any force thatsatisfies the periodicity requirement can be expressed as a sum overharmonics, $$\begin{equation}F(t)=\frac{f_0}{2}+\sum_{n>0} f_n\cos(2n\pi t/\tau)+g_n\sin(2n\pi t/\tau).\label{_auto23} \tag{31}\end{equation}$$ From the previous subsection, one can write down the answer for$x_{pn}(t)$, by substituting $f_n/m$ or $g_n/m$ for $F_0/m$ into Eq.s([16](eq:fastdriven1)) or ([17](eq:fastdriven2)) respectively. Bywriting each factor $2n\pi t/\tau$ as $n\omega t$, with $\omega\equiv2\pi/\tau$, $$\begin{equation}\label{eq:fourierdef1} \tag{32}F(t)=\frac{f_0}{2}+\sum_{n>0}f_n\cos(n\omega t)+g_n\sin(n\omega t).\end{equation}$$ The solutions for $x(t)$ then come from replacing $\omega$ with$n\omega$ for each term in the particular solution in Equations([9](eq:partform)) and ([12](eq:Ddrive)), $$\begin{eqnarray}x_p(t)&=&\frac{f_0}{2k}+\sum_{n>0} \alpha_n\cos(n\omega t-\delta_n)+\beta_n\sin(n\omega t-\delta_n),\\\nonumber\alpha_n&=&\frac{f_n/m}{\sqrt{((n\omega)^2-\omega_0^2)+4\beta^2n^2\omega^2}},\\\nonumber\beta_n&=&\frac{g_n/m}{\sqrt{((n\omega)^2-\omega_0^2)+4\beta^2n^2\omega^2}},\\\nonumber\delta_n&=&\tan^{-1}\left(\frac{2\beta n\omega}{\omega_0^2-n^2\omega^2}\right).\end{eqnarray}$$ Because the forces have been applied for a long time, any non-zerodamping eliminates the homogenous parts of the solution, so one needonly consider the particular solution for each $n$.The problem will considered solved if one can find expressions for thecoefficients $f_n$ and $g_n$, even though the solutions are expressedas an infinite sum. The coefficients can be extracted from thefunction $F(t)$ by $$\begin{eqnarray}\label{eq:fourierdef2} \tag{33}f_n&=&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~F(t)\cos(2n\pi t/\tau),\\\nonumberg_n&=&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~F(t)\sin(2n\pi t/\tau).\end{eqnarray}$$ To check the consistency of these expressions and to verifyEq. ([33](eq:fourierdef2)), one can insert the expansion of $F(t)$ inEq. ([32](eq:fourierdef1)) into the expression for the coefficients inEq. ([33](eq:fourierdef2)) and see whether $$\begin{eqnarray}f_n&=?&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~\left\{\frac{f_0}{2}+\sum_{m>0}f_m\cos(m\omega t)+g_m\sin(m\omega t)\right\}\cos(n\omega t).\end{eqnarray}$$ Immediately, one can throw away all the terms with $g_m$ because theyconvolute an even and an odd function. The term with $f_0/2$disappears because $\cos(n\omega t)$ is equally positive and negativeover the interval and will integrate to zero. For all the terms$f_m\cos(m\omega t)$ appearing in the sum, one can use angle additionformulas to see that $\cos(m\omega t)\cos(n\omegat)=(1/2)(\cos[(m+n)\omega t]+\cos[(m-n)\omega t]$. This will integrateto zero unless $m=n$. In that case the $m=n$ term gives $$\begin{equation}\int_{-\tau/2}^{\tau/2}dt~\cos^2(m\omega t)=\frac{\tau}{2},\label{_auto24} \tag{34}\end{equation}$$ and $$\begin{eqnarray}f_n&=?&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~f_n/2\\\nonumber&=&f_n~\checkmark.\end{eqnarray}$$ The same method can be used to check for the consistency of $g_n$.Consider the driving force: $$\begin{equation}F(t)=At/\tau,~~-\tau/2<t<\tau/2,~~~F(t+\tau)=F(t).\label{_auto25} \tag{35}\end{equation}$$ Find the Fourier coefficients $f_n$ and $g_n$ for all $n$ using Eq. ([33](eq:fourierdef2)).Only the odd coefficients enter by symmetry, i.e. $f_n=0$. One can find $g_n$ integrating by parts, $$\begin{eqnarray}\label{eq:fouriersolution} \tag{36}g_n&=&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2}dt~\sin(n\omega t) \frac{At}{\tau}\\\nonumberu&=&t,~dv=\sin(n\omega t)dt,~v=-\cos(n\omega t)/(n\omega),\\\nonumberg_n&=&\frac{-2A}{n\omega \tau^2}\int_{-\tau/2}^{\tau/2}dt~\cos(n\omega t)+\left.2A\frac{-t\cos(n\omega t)}{n\omega\tau^2}\right|_{-\tau/2}^{\tau/2}.\end{eqnarray}$$ The first term is zero because $\cos(n\omega t)$ will be equallypositive and negative over the interval. Using the fact that$\omega\tau=2\pi$, $$\begin{eqnarray}g_n&=&-\frac{2A}{2n\pi}\cos(n\omega\tau/2)\\\nonumber&=&-\frac{A}{n\pi}\cos(n\pi)\\\nonumber&=&\frac{A}{n\pi}(-1)^{n+1}.\end{eqnarray}$$  Fourier SeriesMore text will come here, chpater 5.7-5.8 of Taylor are discussedduring the lectures. The code here uses the Fourier series discussedin chapter 5.7 for a square wave signal. The equations for thecoefficients are are discussed in Taylor section 5.7, see Example5.4. The code here visualizes the various approximations given byFourier series compared with a square wave with period $T=0.2$, witth$0.1$ and max value $F=2$. We see that when we increase the number ofcomponents in the Fourier series, the Fourier series approximation gets closes and closes to the square wave signal.<jupyter_code>import numpy as np
import math
from scipy import signal
import matplotlib.pyplot as plt

# number of points                                                                                       
n = 500
# start and final times                                                                                  
t0 = 0.0
tn = 1.0
# Period                                                                                                 
T =0.2
# Max value of square signal                                                                             
Fmax= 2.0
# Width of signal                                                                                        
Width = 0.1
t = np.linspace(t0, tn, n, endpoint=False)
SqrSignal = np.zeros(n)
FourierSeriesSignal = np.zeros(n)
SqrSignal = 1.0+signal.square(2*np.pi*5*t+np.pi*Width/T)
a0 = Fmax*Width/T
FourierSeriesSignal = a0
Factor = 2.0*Fmax/np.pi
for i in range(1,500):
    FourierSeriesSignal += Factor/(i)*np.sin(np.pi*i*Width/T)*np.cos(i*t*2*np.pi/T)
plt.plot(t, SqrSignal)
plt.plot(t, FourierSeriesSignal)
plt.ylim(-0.5, 2.5)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15329678470884575, 'fraction_numerical': 0.025514891094977035, 'mean_word_length': 4.41669341894061, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 46, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '11428', 'license': 'CC0-1.0', 'n_tokens_mistral': 11591, 'n_tokens_neox': 11044, 'n_words': 3469, 'path': 'doc/pub/week9/ipynb/.ipynb_checkpoints/week9-checkpoint.ipynb', 'repo_name': 'Shield94/Physics321', 'size': '33773'}