@@ -38,7 +38,7 @@ point for method development, implementation, and analysis.
3838
3939===== A basic model for vibrations =====
4040
41- idx{vibration ODE} idx{oscillations} idx{mechanical vibrations}
41+ idx{vibration ODE} idx{oscillations} idx{mechanical vibrations} idx{angular frequency}
4242
4343The simplest model of a vibrating mechanical system has the following form:
4444
@@ -90,7 +90,7 @@ To formulate a finite difference method for the model
9090problem (ref{vib:ode1}) we follow the ref[four steps explained in Section
9191ref{decay:schemes:FE}][ in cite{Langtangen_decay}]["four steps": "${doc_notes}/sphinx-decay/main_decay.html#the-forward-euler-scheme" explained in cite{Langtangen_decay}].
9292
93- idx{mesh!finite differences} idx{mesh function}
93+ idx{mesh!finite differences} idx{mesh function} idx{discretization of domain}
9494
9595=== Step 1: Discretizing the domain ===
9696
@@ -184,6 +184,8 @@ method, "Verlet integration":
184184that Leapfrog is used for many quite different methods for quite
185185different differential equations!).
186186
187+ idx{Verlet integration} idx{Stoermer's method} idx{Leapfrog method}
188+
187189=== Computing the first step ===
188190
189191We observe that (ref{vib:ode1:step4}) cannot be used for $n=0$ since
@@ -236,6 +238,8 @@ and comparison with the mathematical $\omega$ instead of
236238the full word `omega` as variable name.
237239!ewarning
238240
241+ idx{FD operator notation}
242+
239243=== Operator notation ===
240244
241245We may write the scheme using a compact difference notation
@@ -412,7 +416,7 @@ idx{pytest}
412416idx{nose}
413417idx{verification!hand calculations}
414418idx{unit testing}
415-
419+ idx{zeros}
416420
417421=== Manual calculation ===
418422
@@ -457,7 +461,8 @@ vib_undamped.py::test_convergence_rates PASSED
457461=========================== 2 passed in 0.19 seconds ===...
458462!ec
459463
460- idx{verification!polynomial solutions}
464+ idx{verification!polynomial solutions}
465+ idx{assert}
461466
462467=== Testing very simple polynomial solutions ===
463468
@@ -477,6 +482,7 @@ of the discrete equations. Such results are very useful for debugging
477482and verification. You are strongly encouraged to do this problem now!
478483
479484idx{verification!convergence rates}
485+ idx{error norm}
480486
481487=== Checking convergence rates ===
482488
@@ -498,7 +504,7 @@ In the present problem, computing convergence rates means that we must
498504 $E_i=\sqrt{\Delta t_i\sum_{n=0}^{N_t-1}(u^n-\uex(t_n))^2}$ in each case;
499505 * estimate the convergence rates $r_i$ based on two consecutive
500506 experiments $(\Delta t_{i-1}, E_{i-1})$ and $(\Delta t_{i}, E_{i})$,
501- assuming $E_i=C(\Delta t_i)^{r}$ and $E_{i-1}=C(\Delta t_{i-1})^{r}$.
507+ assuming $E_i=C(\Delta t_i)^{r}$ and $E_{i-1}=C(\Delta t_{i-1})^{r}$, where $C$ is a constant .
502508 From these equations it follows that
503509 $r = \ln (E_{i-1}/E_i)/\ln (\Delta t_{i-1}/\Delta t_i)$. Since this $r$
504510 will vary with $i$, we equip it with an index and call it $r_{i-1}$,
@@ -540,8 +546,9 @@ error model $E_i = C(\Delta t_i)^r$ is valid.
540546
541547We can now construct a proper test function that computes convergence rates
542548and checks that the final (and usually the best) estimate is sufficiently
543- close to 2. Here, a rough tolerance of 0.1 is enough. This unit test
544- goes like
549+ close to 2. Here, a rough tolerance of 0.1 is enough. Later, we will argue
550+ for an improvement by adjusting omega and include also that case in our test
551+ function here. The unit test goes like
545552
546553@@@CODE src-vib/vib_undamped.py fromto: def test_convergence@def plot_convergence_rates
547554The complete code appears in the file `vib_undamped.py`.
@@ -645,6 +652,8 @@ FIGURE: [fig-vib/vib_freq_err1, width=800 frac=1.0] Effect of halving the time s
645652
646653===== Using a moving plot window =====
647654
655+ idx{SciTools}
656+
648657In vibration problems it is often of interest to investigate the system's
649658behavior over long time intervals. Errors in the angular frequency accumulate
650659and become more visible as time grows. We can investigate long
@@ -796,6 +805,8 @@ Terminal> mv vib.html vib_dt0.1/index.html
796805
797806=== Making animated GIF files ===
798807
808+ idx{ImageMagic} idx{Bokeh}
809+
799810The `convert` program from the ImageMagick software suite can be
800811used to produce animated GIF files from a set of PNG files:
801812
@@ -849,7 +860,7 @@ FIGURE: [fig-vib/bokeh_gridplot_interactive, width=800 frac=1] Interactive Bokeh
849860
850861A function for creating a Bokeh plot, given a list of `u` arrays
851862and corresponding `t` arrays, is implemented below.
852- The code combines data fro different simulations, described
863+ The code combines data from different simulations, described
853864compactly in a list of strings `legends`.
854865
855866@@@CODE src-vib/vib_undamped.py fromto: def bokeh_plot@def demo_bokeh
@@ -1109,7 +1120,9 @@ label{vib:ode1:tildeomega}
11091120===== The error in the numerical frequency =====
11101121label{vib:ode1:analysis:numfreq}
11111122
1112- The first observation of (ref{vib:ode1:tildeomega}) tells that there
1123+ idx{sympy}
1124+
1125+ The first observation following (ref{vib:ode1:tildeomega}) tells that there
11131126is a phase error since the numerical frequency $\tilde\omega$ never
11141127equals the exact frequency $\omega$. But how good is the approximation
11151128(ref{vib:ode1:tildeomega})? That is, what is the error $\omega -
@@ -1171,7 +1184,7 @@ FIGURE: [fig-vib/discrete_freq, width=400] Exact discrete frequency and its seco
11711184
11721185===== Empirical convergence rates and adjusted $\omega$ =====
11731186
1174- The expression (ref{vib:ode1:tildeomega:series}) suggest that
1187+ The expression (ref{vib:ode1:tildeomega:series}) suggests that
11751188adjusting omega to
11761189
11771190!bt
@@ -1199,6 +1212,8 @@ impact on $\omega$, and therefore we cannot use the trick.
11991212===== Exact discrete solution =====
12001213label{vib:ode1:analysis:sol}
12011214
1215+ idx{error mesh function}
1216+
12021217Perhaps more important than the $\tilde\omega = \omega + {\cal O}(\Delta t^2)$
12031218result found above is the fact that we have an exact discrete solution of
12041219the problem:
@@ -1240,8 +1255,8 @@ by doing Exercise ref{vib:exer:discrete:omega}.
12401255label{vib:ode1:analysis:conv}
12411256
12421257
1243- We can use (ref{vib:ode1:tildeomega:series}), (ref{vib:ode1:en}), or
1244- (ref{vib:ode1:en2}) to show *convergence* of the numerical scheme,
1258+ We can use (ref{vib:ode1:tildeomega:series}) and (ref{vib:ode1:en}), or
1259+ (ref{vib:ode1:en2}), to show *convergence* of the numerical scheme,
12451260i.e., $e^n\rightarrow 0$ as $\Delta t\rightarrow 0$, which implies
12461261that the numerical solution approaches the exact solution as $\Delta
12471262t$ approaches to zero. We have that
@@ -1269,7 +1284,7 @@ It then follows from the expression(s) for $e^n$ that $e^n\rightarrow 0$.
12691284
12701285===== The global error =====
12711286
1272- idx{error!global}
1287+ idx{error!global} idx{error norm} idx{norm}
12731288
12741289To achieve more analytical insight into the nature of the global
12751290error, we can Taylor expand the error mesh function
@@ -1458,6 +1473,8 @@ analysis we can draw three important conclusions:
14581473======= Alternative schemes based on 1st-order equations =======
14591474label{vib:model2x2}
14601475
1476+ idx{1st-order ODE} idx{2nd-order ODE}
1477+
14611478A standard technique for solving second-order ODEs is to rewrite them
14621479as a system of first-order ODEs and then choose a solution strategy
14631480from the vast collection of methods for first-order ODE systems.
@@ -1691,6 +1708,8 @@ terms at the same time point $t_n$.
16911708===== Comparison of schemes =====
16921709label{vib:model2x2:compare}
16931710
1711+ idx{Odespy}
1712+
16941713We can easily compare methods like the ones above (and many more!)
16951714with the aid of the
16961715"Odespy": "https://github.com/hplgit/odespy" package. Below is
@@ -1932,6 +1951,9 @@ dive into this subject.
19321951===== Derivation of the energy expression =====
19331952label{vib:model1:energy:expr}
19341953
1954+ idx{kinetic energy} idx{potential energy} idx{spring constant} idx{stiffness}
1955+ idx{Newton's 2nd law}
1956+
19351957We start out with multiplying
19361958
19371959!bt
@@ -2039,8 +2061,8 @@ is constant:
20392061
20402062=== Growth of energy in the Forward Euler scheme ===
20412063
2042- The energy at time level $n+1$ in the Forward Euler scheme can easily
2043- be shown to increase:
2064+ It is easy to show that the energy in the Forward Euler scheme increases
2065+ when stepping from time level $n$ to $n+1$.
20442066
20452067!bt
20462068\begin{align*}
@@ -2134,6 +2156,8 @@ is subtracted.
21342156===== An error measure based on energy =====
21352157label{vib:model1:energy:measure}
21362158
2159+ idx{error norm}
2160+
21372161The constant energy is well expressed by its initial value $E(0)$, so that
21382162the error in mechanical energy can be computed as a mesh function by
21392163
@@ -2309,12 +2333,14 @@ the second-order equation $u^{\prime\prime}+\omega^2u=0$. However,
23092333there is a similarly successful scheme available for the first-order
23102334system $u^{\prime}=v$, $v^{\prime}=-\omega^2u$, to be presented below.
23112335The ideas of the scheme and their further developments have become
2312- very popular in particle and rigid body dynamics and hence widely
2336+ very popular in particle and rigid body dynamics and hence are widely
23132337used by physicists.
23142338
23152339
2316- idx{forward-backward Euler-Cromer scheme}
2340+ idx{forward-backward scheme}
23172341idx{Euler-Cromer scheme}
2342+ idx{semi-implicit Euler}
2343+ idx{semi-explicit Euler}
23182344
23192345===== Forward-backward discretization =====
23202346
@@ -2406,12 +2432,6 @@ Newton-St\"{o}rmer-Verlet,
24062432Newton-Stoermer-Verlet,
24072433# #endif
24082434and Euler-Cromer. We shall stick to the latter name.
2409- Since both time
2410- discretizations are based on first-order difference approximation, one
2411- may think that the scheme is only of first-order, but this is not
2412- true: the use of a forward and then a backward difference make errors
2413- cancel so that the overall error in the scheme is $\Oof{\Delta
2414- t^2}$. This is explained below. [hpl: This is not correct!]
24152435
24162436How does the Euler-Cromer method preserve the total energy?
24172437We may run the example from Section ref{vib:model1:energy:measure}:
@@ -2474,7 +2494,7 @@ method.
24742494
24752495Although the Euler-Cromer scheme and the method (ref{vib:ode1:step4})
24762496are equivalent, there could be differences in the way they handle the
2477- initial conditions. Let is look into this topic. The initial
2497+ initial conditions. Let us look into this topic. The initial
24782498condition $u^{\prime}=0$ means $u^{\prime}=v=0$. From
24792499(ref{vib:model2x2:EulerCromer:veq1b}) we get
24802500
@@ -2586,15 +2606,17 @@ u = u[:,1] # Extract displacement
25862606
25872607=== Convergence rates ===
25882608
2609+ idx{verification! convergence rates}
2610+
25892611We may use the `convergence_rates` function in the file
25902612`vib_undamped.py` to investigate the convergence rate of the
25912613Euler-Cromer method, see the `convergence_rate` function in the file
25922614`vib_undamped_EulerCromer.py`. Since we could eliminate $v$ to get a
25932615scheme for $u$ that is equivalent to the finite difference method for
25942616the second-order equation in $u$, we would expect the convergence
2595- rates to be the same, i.e., $\mathcal{O}(\Delta t^2) $. However,
2617+ rates to be the same, i.e., $r = 2 $. However,
25962618measuring the convergence rate of $u$ in the Euler-Cromer scheme shows
2597- that it is $\mathcal{O}(\Delta t)$ ! Adjusting the initial condition
2619+ that $r = 1$ only ! Adjusting the initial condition
25982620does not change the rate. Adjusting $\omega$, as outlined in Section
25992621ref{vib:ode1:analysis:numfreq}, gives a 4th-order method there, while
26002622there is no increase in the measured rate in the Euler-Cromer
@@ -2622,7 +2644,7 @@ idx{Stoermer-Verlet algorithm}
26222644Another very popular algorithm for vibration problems, especially
26232645for long time simulations, is the
26242646# #if FORMAT in ("pdflatex", "latex")
2625- St\"{o}mer -Verlet
2647+ St\"{o}rmer -Verlet
26262648# #else
26272649Stoermer-Verlet
26282650# #endif
@@ -2643,14 +2665,14 @@ With mathematics,
26432665\begin{align*}
26442666\frac{v^{n+\half}-v^n}{\half\Delta t} &= -\omega^2 u^n,\\
26452667\frac{u^{n+\half}-u^n}{\half\Delta t} &= v^{n+\half},\\
2646- \frac{u^{n+1}-u^{n- \half}}{\half\Delta t} &= v^{n+\half},\\
2668+ \frac{u^{n+1}-u^{n+ \half}}{\half\Delta t} &= v^{n+\half},\\
26472669\frac{v^{n+1}-v^{n+\half}}{\half\Delta t} &= -\omega^2 u^{n+1}\tp
26482670\end{align*}
26492671!et
26502672The two steps in the middle can be combined to
26512673
26522674!bt
2653- \[ \frac{u^{n+1}-u^{n-1 }}{\Delta t} = v^{n+\half},\]
2675+ \[ \frac{u^{n+1}-u^{n}}{\Delta t} = v^{n+\half},\]
26542676!et
26552677and consequently
26562678
@@ -2684,7 +2706,7 @@ St\"{o}mer-Verlet,
26842706# #else
26852707Stoermer-Verlet,
26862708# #endif
2687- and using centered differences for the $2times 2$ system on a staggered
2709+ and using centered differences for the $2\times 2$ system on a staggered
26882710mesh) all end up with the same equations! The main difference is that
26892711the traditional Euler-Cromer displays first-order convergence in $\Delta t$
26902712(due to less symmetry in the way $u$ and $v$ are treated)
@@ -2700,7 +2722,7 @@ of explaining it, we must give in too. At least we have shown that
27002722the popular Stoermer-Verlet method is nothing but our staggered mesh very
27012723natural scheme.]
27022724
2703- The most numerical stable scheme, with respect to accumulation of
2725+ The most numerically stable scheme, with respect to accumulation of
27042726rounding errors, is
27052727(ref{vib:model2x2:StormerVerlet:eqv})-(ref{vib:model2x2:StormerVerlet:equ}).
27062728It has, according to cite{Hairer_Wanner_Norsett_bookI}, better
0 commit comments