From eacd6f7bc1a4a048e1352b740dd3354e2a035106 Mon Sep 17 00:00:00 2001 From: Prefetch Date: Fri, 11 Mar 2022 21:15:23 +0100 Subject: Expand knowledge base --- .../know/concept/archimedes-principle/index.pdc | 11 +- content/know/concept/metacentric-height/index.pdc | 185 ++++++++++++++ content/know/concept/metacentric-height/sketch.png | Bin 0 -> 171361 bytes content/know/concept/runge-kutta-method/index.pdc | 267 +++++++++++++++++++++ 4 files changed, 457 insertions(+), 6 deletions(-) create mode 100644 content/know/concept/metacentric-height/index.pdc create mode 100644 content/know/concept/metacentric-height/sketch.png create mode 100644 content/know/concept/runge-kutta-method/index.pdc (limited to 'content') diff --git a/content/know/concept/archimedes-principle/index.pdc b/content/know/concept/archimedes-principle/index.pdc index 3a063ec..fb91b67 100644 --- a/content/know/concept/archimedes-principle/index.pdc +++ b/content/know/concept/archimedes-principle/index.pdc @@ -39,7 +39,7 @@ $$\begin{aligned} Where $\va{g}$ is the gravitational field, and $\rho_\mathrm{b}$ is the density of the body. Meanwhile, the pressure $p$ of the surrounding fluid exerts a force -on the surface $S$ of $V$: +on the entire surface $S$ of $V$: $$\begin{aligned} \va{F}_p @@ -75,18 +75,17 @@ and zero on the "non-submerged" side, we find: $$\begin{aligned} 0 - = \mathrm{g} (\rho_\mathrm{b} - \rho_\mathrm{f}) V = \mathrm{g} (m_\mathrm{b} - m_\mathrm{f}) \end{aligned}$$ -In other words, the mass $m_\mathrm{b}$ of the submerged portion $V$ of the body, +In other words, the mass $m_\mathrm{b}$ of the entire body is equal to the mass $m_\mathrm{f}$ of the fluid it displaces. This is the best-known version of Archimedes' principle. -Note that if $\rho_\mathrm{b} > \rho_\mathrm{f}$, then, +Note that if $\rho_\mathrm{b} > \rho_\mathrm{f}$, +then the displaced mass $m_\mathrm{f} < m_\mathrm{b}$ even if the entire body is submerged, -the displaced mass $m_\mathrm{f} < m_\mathrm{b}$, -and the object will continue to sink. +and the object will therefore continue to sink. diff --git a/content/know/concept/metacentric-height/index.pdc b/content/know/concept/metacentric-height/index.pdc new file mode 100644 index 0000000..1fc6aca --- /dev/null +++ b/content/know/concept/metacentric-height/index.pdc @@ -0,0 +1,185 @@ +--- +title: "Metacentric height" +firstLetter: "M" +publishDate: 2022-03-11 +categories: +- Physics +- Fluid mechanics + +date: 2021-05-08T19:03:36+02:00 +draft: false +markup: pandoc +--- + +# Metacentric height + +Consider an object with center of mass $G$, +floating in a large body of liquid whose surface is flat at $z = 0$. +For our purposes, it is easiest to use a coordinate system +whose origin is at the area centroid +of the object's cross-section through the liquid's surface, namely: + +$$\begin{aligned} + (x_0, y_0) + \equiv \frac{1}{A_{wl}} \iint_{wl} (x, y) \dd{A} +\end{aligned}$$ + +Where $A_{wl}$ is the cross-sectional area +enclosed by the "waterline" around the "boat". +Note that the boat's center of mass $G$ +does not coincide with the origin in general, +as is illustrated in the following sketch +of our choice of coordinate system: + + + + + +Here, $B$ is the **center of buoyancy**, equal to +the center of mass of the volume of water displaced by the boat +as per [Archimedes' principle](/know/concept/archimedes-principle/). +At equilibrium, the forces of buoyancy $\vb{F}_B$ and gravity $\vb{F}_G$ +have equal magnitudes in opposite directions, +and $B$ is directly above or below $G$, +or in other words, $x_B = x_G$ and $y_B = y_G$, +which are calculated as follows: + +$$\begin{aligned} + (x_G, y_G, z_G) + &\equiv \frac{1}{V_{boat}} \iiint_{boat} (x, y, z) \dd{V} + \\ + (x_B, y_B, z_B) + &\equiv \frac{1}{V_{disp}} \iiint_{disp} (x, y, z) \dd{V} +\end{aligned}$$ + +Where $V_{boat}$ is the volume of the whole boat, +and $V_{disp}$ is the volume of liquid it displaces. + +Whether a given equilibrium is *stable* is more complicated. +Suppose the ship is tilted by a small angle $\theta$ around the $x$-axis, +in which case the old waterline, previously in the $z = 0$ plane, +gets shifted to a new plane, namely: + +$$\begin{aligned} + z + = \sin\!(\theta) \: y + \approx \theta y +\end{aligned}$$ + +Then $V_{disp}$ changes by $\Delta V_{disp}$, which is estimated below. +If a point of the old waterline is raised by $z$, +then the displaced liquid underneath it is reduced proportionally, +hence the sign: + +$$\begin{aligned} + \Delta V_{disp} + \approx - \iint_{wl} z \dd{A} + \approx - \theta \iint_{wl} y \dd{A} + = 0 +\end{aligned}$$ + +So $V_{disp}$ is unchanged, at least to first order in $\theta$. +However, the *shape* of the displaced volume may have changed significantly. +Therefore, the shift of the position of the buoyancy center from $B$ to $B'$ +involves a correction $\Delta y_B$ in addition to the rotation by $\theta$: + +$$\begin{aligned} + y_B' + = y_B - \theta z_B + \Delta y_B +\end{aligned}$$ + +We find $\Delta y_B$ by calculating the virtual buoyancy center of the shape difference: +on the side of the boat that has been lifted by the rotation, +the center of buoyancy is "pushed" away due to the reduced displacement there, +and vice versa on the other side. Consequently: + +$$\begin{aligned} + \Delta y_B + = - \frac{1}{V_{disp}} \iint_{wl} y z \dd{A} + \approx - \frac{\theta}{V_{disp}} \iint_{wl} y^2 \dd{A} + = - \frac{\theta I}{V_{disp}} +\end{aligned}$$ + +Where we have defined the so-called **area moment** $I$ of the waterline as follows: + +$$\begin{aligned} + \boxed{ + I + \equiv \iint_{wl} y^2 \dd{A} + } +\end{aligned}$$ + +Now that we have an expression for $\Delta y_B$, +the new center's position $y_B'$ is found to be: + +$$\begin{aligned} + y_B' + = y_B - \theta \Big( z_B + \frac{I}{V_{disp}} \Big) + \approx y_B - \sin\!(\theta) \: \Big( z_B + \frac{I}{V_{disp}} \Big) +\end{aligned}$$ + +This looks like a rotation by $\theta$ around a so-called **metacenter** $M$, +with a height $z_M$ known as the **metacentric height**, defined as: + +$$\begin{aligned} + \boxed{ + z_M + \equiv z_B + \frac{I}{V_{disp}} + } +\end{aligned}$$ + +Meanwhile, the position of $M$ is defined such that it lies +on the line between the old centers $G$ and $B$. +Our calculation of $y_B'$ has shown that the new $B'$ always lies below $M$. + +After the rotation, the boat is not in equilibrium anymore, +because the new $G'$ is not directly above or below $B'$. +The force of gravity then causes a torque $\vb{T}$ given by: + +$$\begin{aligned} + \vb{T} + = (\vb{r}_G' - \vb{r}_B') \cross m \vb{g} +\end{aligned}$$ + +Where $\vb{g}$ points downwards. +Since the rotation was around the $x$-axis, +we are only interested in the $x$-component $T_x$, which becomes: + +$$\begin{aligned} + T_x + = - (y_G' - y_B') m \mathrm{g} + = - \big((y_G - \theta z_G) - (y_B - \theta z_M)\big) m \mathrm{g} +\end{aligned}$$ + +With $y_G' = y_G - \theta z_G$ being a simple rotation of $G$. +At the initial equilibrium $y_G = y_B$, so: + +$$\begin{aligned} + T_x + = \theta (z_G - z_M) m \mathrm{g} +\end{aligned}$$ + +If $z_M < z_G$, then $T_x$ has the same sign as $\theta$, +so $\vb{T}$ further destabilizes the boat. +But if $z_M > z_G$, then $\vb{T}$ counteracts the rotation, +and the boat returns to the original equilibrium, +leading us to the following stability condition: + +$$\begin{aligned} + \boxed{ + z_M > z_G + } +\end{aligned}$$ + +In other words, for a given boat design (or general shape) +$z_G$ and $z_M$ can be calculated, +and as long as they satisfy the above inequality, +it will float stably in water (or any other fluid, +although the buoyancy depends significantly on the density). + + + +## References +1. B. Lautrup, + *Physics of continuous matter: exotic and everyday phenomena in the macroscopic world*, 2nd edition, + CRC Press. diff --git a/content/know/concept/metacentric-height/sketch.png b/content/know/concept/metacentric-height/sketch.png new file mode 100644 index 0000000..f33fe36 Binary files /dev/null and b/content/know/concept/metacentric-height/sketch.png differ diff --git a/content/know/concept/runge-kutta-method/index.pdc b/content/know/concept/runge-kutta-method/index.pdc new file mode 100644 index 0000000..ac2eabf --- /dev/null +++ b/content/know/concept/runge-kutta-method/index.pdc @@ -0,0 +1,267 @@ +--- +title: "Runge-Kutta method" +firstLetter: "R" +publishDate: 2022-03-10 +categories: +- Mathematics +- Numerical methods + +date: 2022-03-07T14:10:18+01:00 +draft: false +markup: pandoc +--- + +# Runge-Kutta method + +A **Runge-Kutta method** (RKM) is a popular approach +to numerically solving systems of ordinary differential equations. +Let $\vb{x}(t)$ be the vector we want to find, +governed by $\vb{f}(t, \vb{x})$: + +$$\begin{aligned} + \vb{x}'(t) + = \vb{f}\big(t, \vb{x}(t)\big) +\end{aligned}$$ + +Like in all numerical methods, the $t$-axis is split into discrete steps. +If a step has size $h$, then as long as $h$ is small enough, +we can make the following approximation: + +$$\begin{aligned} + \vb{x}'(t) + a h \vb{x}''(t) + &\approx \vb{x}'(t \!+\! a h) + \\ + &\approx \vb{f}\big(t \!+\! a h,\, \vb{x}(t \!+\! a h)\big) + \\ + &\approx \vb{f}\big(t \!+\! a h,\, \vb{x}(t) \!+\! a h \vb{x}'(t) \big) +\end{aligned}$$ + +For sufficiently small $h$, +higher-order derivates can also be included, +albeit still at $t \!+\! a h$: + +$$\begin{aligned} + \vb{x}'(t) + a h \vb{x}''(t) + b h^2 \vb{x}'''(t) + &\approx \vb{f}\big(t \!+\! a h,\, \vb{x}(t) \!+\! a h \vb{x}'(t) \!+\! b h^2 \vb{x}''(t) \big) +\end{aligned}$$ + +Although these approximations might seem innocent, +they actually make it quite complicated to determine the error order of a given RKM. + +Now, consider a Taylor expansion around the current $t$, +truncated at a chosen order $n$: + +$$\begin{aligned} + \vb{x}(t \!+\! h) + &= \vb{x}(t) + h \vb{x}'(t) + \frac{h^2}{2} \vb{x}''(t) + \frac{h^3}{6} \vb{x}'''(t) + \:...\, + \frac{h^n}{n!} \vb{x}^{(n)}(t) + \\ + &= \vb{x}(t) + h \bigg[ \vb{x}'(t) + \frac{h}{2} \vb{x}''(t) + \frac{h^2}{6} \vb{x}'''(t) + \:...\, + \frac{h^{n-1}}{n!} \vb{x}^{(n)}(t) \bigg] +\end{aligned}$$ + +We are free to split the terms as follows, +choosing real factors $\omega_{mj}$ subject to $\sum_{j} \omega_{mj} = 1$: + +$$\begin{aligned} + \vb{x}(t \!+\! h) + &= \vb{x} + h \bigg[ \sum_{j = 1}^{N_1} \omega_{1j} \, \vb{x}' + + \frac{h}{2} \sum_{j = 1}^{N_2} \omega_{2j} \, \vb{x}'' + + \:...\, + \frac{h^{n-1}}{n!} \sum_{j = 1}^{N_n} \omega_{nj} \, \vb{x}^{(n)} \bigg] +\end{aligned}$$ + +Where the integers $N_1,...,N_n$ are also free to choose, +but for reasons that will become clear later, +the most general choice for an RKM is $N_1 = n$, $N_n = 1$, and: + +$$\begin{aligned} + N_{n-1} + = N_n \!+\! 2 + ,\quad + \cdots + ,\quad + N_{n-m} + = N_{n-m+1} \!+\! m \!+\! 1 + ,\quad + \cdots + ,\quad + N_{2} + = N_3 \!+\! n \!-\! 1 +\end{aligned}$$ + +In other words, $N_{n-m}$ is the $m$th triangular number. +This is not so important, +since this is not a practical way to describe RKMs, +but it is helpful to understand how they work. + + +## Example derivation + +For example, let us truncate at $n = 3$, +such that $N_1 = 3$, $N_2 = 3$ and $N_3 = 1$. +The following derivation is very general, +except it requires all $\alpha_j \neq 0$. +Renaming $\omega_{mj}$, we start from: + +$$\begin{aligned} + \vb{x}(t \!+\! h) + &= \vb{x} + h \bigg[ (\alpha_1 + \alpha_2 + \alpha_3) \, \vb{x}' + + \frac{h}{2} (\beta_2 + \beta_{31} + \beta_{32}) \, \vb{x}'' + + \frac{h^2}{6} \gamma_3 \, \vb{x}''' \bigg] + \\ + &= \vb{x} + h \bigg[ \alpha_1 \vb{x}' + + \Big( \alpha_2 \vb{x}' + \frac{h}{2} \beta_2 \vb{x}'' \Big) + + \Big( \alpha_3 \vb{x}' + \frac{h}{2} (\beta_{31} + \beta_{32}) \vb{x}'' + \frac{h^2}{6} \gamma_3 \vb{x}''' \Big) \bigg] +\end{aligned}$$ + +As discussed earlier, the parenthesized expressions +can be approximately rewritten with $\vb{f}$: + +$$\begin{aligned} + \vb{x}(t \!+\! h) + = \vb{x} + h &\bigg[ \alpha_1 \vb{f}(t, \vb{x}) + + \alpha_2 \vb{f}\Big( t \!+\! \frac{h \beta_2}{2 \alpha_2}, \; + \vb{x} \!+\! \frac{h \beta_2}{2 \alpha_2} \vb{x}' \Big) + \\ + & + \alpha_3 \vb{f}\Big( t \!+\! \frac{h (\beta_{31} \!\!+\!\! \beta_{32})}{2 \alpha_3}, \; + \vb{x} \!+\! \frac{h \beta_{31}}{2 \alpha_3} \vb{x}' \!+\! \frac{h \beta_{32}}{2 \alpha_3} \vb{x}' + \!+\! \frac{h^2 \gamma_3}{6 \alpha_3} \vb{x}'' \Big) \bigg] + \\ + = \vb{x} + h &\bigg[ \alpha_1 \vb{k}_1 + + \alpha_2 \vb{f}\Big( t \!+\! \frac{h \beta_2}{2 \alpha_2}, \; + \vb{x} \!+\! \frac{h \beta_2}{2 \alpha_2} \vb{k}_1 \!\Big) + \\ + & + \alpha_3 \vb{f}\Big( t \!+\! \frac{h (\beta_{31} \!\!+\!\! \beta_{32})}{2 \alpha_3}, \; + \vb{x} \!+\! \frac{h \beta_{31}}{2 \alpha_3} \vb{k}_1 \!+\! \frac{h \beta_{32}}{2 \alpha_3} + \vb{f}\Big( t \!+\! \frac{h \gamma_3}{3 \beta_{32}}, \; + \vb{x} \!+\! \frac{h \gamma_3}{3 \beta_{32}} \vb{k}_1 \!\Big) \!\Big) \bigg] +\end{aligned}$$ + +Here, we can see an opportunity to save some computational time +by reusing an evaluation of $\vb{f}$. +Technically, this is optional, but it would be madness not to, +so we choose: + +$$\begin{aligned} + \frac{\beta_2}{2 \alpha_2} + = \frac{\gamma_3}{3 \beta_{32}} +\end{aligned}$$ + +Such that the next step of $\vb{x}$'s numerical solution is as follows, +recalling that $\sum_{j} \alpha_j = 1$: + +$$\begin{aligned} + \boxed{ + \vb{x}(t \!+\! h) + = \vb{x}(t) + h \Big( \alpha_1 \vb{k}_1 + \alpha_2 \vb{k}_2 + \alpha_3 \vb{k}_3 \Big) + } +\end{aligned}$$ + +Where $\vb{k}_1$, $\vb{k}_2$ and $\vb{k}_3$ are different estimates +of the average slope $\vb{x}'$ between $t$ and $t \!+\! h$, +whose weighted average is used to make the $t$-step. +They are given by: + +$$\begin{aligned} + \boxed{ + \begin{aligned} + \vb{k}_1 + &\equiv \vb{f}(t, \vb{x}) + \\ + \vb{k}_2 + &\equiv \vb{f}\bigg( t + \frac{h \beta_2}{2 \alpha_2}, \; + \vb{x} + \frac{h \beta_2}{2 \alpha_2} \vb{k}_1 \bigg) + \\ + \vb{k}_3 + &\equiv \vb{f}\bigg( t + \frac{h (\beta_{31} \!\!+\!\! \beta_{32})}{2 \alpha_3}, \; + \vb{x} + \frac{h \beta_{31}}{2 \alpha_3} \vb{k}_1 + \frac{h \beta_{32}}{2 \alpha_3} \vb{k}_2 \bigg) + \end{aligned} + } +\end{aligned}$$ + +Despite the contraints on $\alpha_j$ and $\beta_j$, +there is an enormous freedom of choice here, +all leading to valid RKMs, although not necessarily good ones. + + +## General form + +A more practical description goes as follows: +in an $s$-stage RKM, a weighted average is taken +of up to $s$ slope estimates $\vb{k}_j$ with weights $b_j$. +Let $\sum_{j} b_j = 1$, then: + +$$\begin{aligned} + \boxed{ + \vb{x}(t \!+\! h) + = \vb{x}(t) + h \sum_{j = 1}^{s} b_j \vb{k}_j + } +\end{aligned}$$ + +Where the estimates $\vb{k}_1, ..., \vb{k}_s$ +depend on each other, and are calculated one by one as: + +$$\begin{aligned} + \boxed{ + \vb{k}_m + = \vb{f}\bigg( t + h c_m,\; \vb{x} + h \sum_{j = 1}^{m - 1} a_{mj} \vb{k}_j \bigg) + } +\end{aligned}$$ + +With $c_1 = 1$ and $\sum_{j = 1} a_{mj} = c_m$. +Writing this out for the first few $m$, the pattern is clear: + +$$\begin{aligned} + \vb{k}_1 + &= \vb{f}(t, \vb{x}) + \\ + \vb{k}_2 + &= \vb{f}\big( t + h c_2,\; \vb{x} + h a_{21} \vb{k}_1 \big) + \\ + \vb{k}_3 + &= \vb{f}\big( t + h c_3,\; \vb{x} + h (a_{31} \vb{k}_1 + a_{32} \vb{k}_2) \big) + \\ + \vb{k}_4 + &= \:... +\end{aligned}$$ + +The coefficients of a given RKM are usually +compactly represented in a **Butcher tableau**: + +$$\begin{aligned} + \begin{array}{c|ccc} + 0 \\ + c_2 & a_{21} \\ + c_3 & a_{31} & a_{32} \\ + \vdots & \vdots & \vdots & \ddots \\ + c_s & a_{s1} & a_{s2} & \cdots & a_{s,s-1} \\ + \hline + & b_1 & b_2 & \cdots & b_{s-1} & b_s + \end{array} +\end{aligned}$$ + +Each RKM has an **order** $p$, +such that the global truncation error is $\mathcal{O}(h^p)$, +i.e. the accumulated difference between the numerical +and the exact solutions is proportional to $h^p$. + +The surprise is that $p$ need not be equal to the Taylor expansion order $n$, +nor the stage count $s$. +Typically, $s = n$ for computational efficiency, but $s \ge n$ is possible in theory. + +The order $p$ of a given RKM is determined by +a complicated set of equations on the coefficients, +and the lowest possible $s$ for a desired $p$ +is in fact only partially known. +For $p \le 4$ the bound is $s \ge p$, +whereas for $p \ge 5$ the only proven bound is $s \ge p \!+\! 1$, +but for $p \ge 7$ no such efficient methods have been found so far. + +If you need an RKM with a certain order, look it up. +There exist many efficient methods for $p \le 4$ where $s = p$, +and although less popular, higher $p$ are also available. + + + +## References +1. J.C. Butcher, + *Numerical methods for ordinary differential equations*, 3rd edition, + Wiley. -- cgit v1.2.3