Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: vib_gen #86

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

WIP: vib_gen #86

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6109a41
Start formatting
jakubbober Sep 22, 2020
d09ed38
Start formatting
jakubbober Sep 22, 2020
6305099
Finish formatting
jakubbober Sep 22, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 44 additions & 34 deletions fdm-devito-notebooks/01_vib/vib_gen.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,19 @@
{
"cells": [
{
"cell_type": "raw",
"cell_type": "markdown",
"metadata": {},
"source": [
"# Generalization: damping, nonlinearities, and excitation\n",
"<div id=\"vib:gen\"></div>\n",
"\n",
"We shall now generalize the simple model problem from Section [Finite difference discretization](vib_undamped.ipynb#vib:model1) to\n",
"include a possibly nonlinear damping term $f(u')$, a possibly nonlinear\n",
"spring (or restoring) force $s(u)$, and some external excitation $F(t)$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
Expand Down Expand Up @@ -239,7 +251,7 @@
"metadata": {},
"source": [
"for some quantity $w$ depending on time. The error in the geometric mean\n",
"approximation is $\\Oof{\\Delta t^2}$, the same as in the\n",
"approximation is $\\mathcal{O}(\\Delta t^2)$, the same as in the\n",
"approximation $u^{\\prime\\prime}\\approx D_tD_tu$. With $w=u^{\\prime}$ it follows\n",
"that"
]
Expand Down Expand Up @@ -464,7 +476,7 @@
"result in exactly the same scheme as in\n",
"([11](#vib:ode2:step3b:quad)) where we\n",
"used a geometric mean and centered differences and committed errors\n",
"of size $\\Oof{\\Delta t^2}$. Therefore, the forward-backward\n",
"of size $\\mathcal{O}(\\Delta t^2)$. Therefore, the forward-backward\n",
"differences in ([15](#vib:ode2:nonlin:fbdiff))\n",
"act in a symmetric way and actually produce a second-order\n",
"accurate discretization of the quadratic damping term.\n",
Expand All @@ -476,7 +488,7 @@
"\n",
"The algorithm arising from the methods in the sections [A centered scheme for linear damping](#vib:ode2:fdm:flin)\n",
"and [A centered scheme for quadratic damping](#vib:ode2:fdm:fquad) is very similar to the undamped case in\n",
"the section [vib:ode1:fdm](#vib:ode1:fdm). The difference is\n",
"the section [A centered finite difference scheme](vib_undamped.ipynb#vib:ode1:fdm). The difference is\n",
"basically a question of different formulas for $u^1$ and\n",
"$u^{n+1}$. This is actually quite remarkable. The equation\n",
"([1](#vib:ode2)) is normally impossible to solve by pen and paper, but\n",
Expand Down Expand Up @@ -546,15 +558,15 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"The complete code resides in the file [`vib.py`](${src_vib}/vib.py).\n",
"The complete code resides in the file [`vib.py`](https://github.com/devitocodes/devito_book/blob/master/fdm-devito-notebooks/01_vib/src-vib/vib.py).\n",
"\n",
"## Verification\n",
"<div id=\"vib:ode2:verify\"></div>\n",
"\n",
"### Constant solution\n",
"\n",
"For debugging and initial verification, a constant solution is often\n",
"very useful. We choose $\\uex(t)=I$, which implies $V=0$.\n",
"very useful. We choose $u_\\text{e}(t)=I$, which implies $V=0$.\n",
"Inserted in the ODE, we get\n",
"$F(t)=s(I)$ for any choice of $f$. Since the discrete derivative\n",
"of a constant vanishes (in particular, $[D_{2t}I]^n=0$,\n",
Expand All @@ -565,9 +577,9 @@
"\n",
"### Linear solution\n",
"\n",
"Now we choose a linear solution: $\\uex = ct + d$. The initial condition\n",
"Now we choose a linear solution: $u_\\text{e} = ct + d$. The initial condition\n",
"$u(0)=I$ implies $d=I$, and $u^{\\prime}(0)=V$ forces $c$ to be $V$.\n",
"Inserting $\\uex=Vt+I$ in the ODE with linear damping results in"
"Inserting $u_\\text{e}=Vt+I$ in the ODE with linear damping results in"
]
},
{
Expand Down Expand Up @@ -600,21 +612,21 @@
"metadata": {},
"source": [
"Since the finite difference approximations used to compute $u^{\\prime}$ all\n",
"are exact for a linear function, it turns out that the linear $\\uex$\n",
"are exact for a linear function, it turns out that the linear $u_\\text{e}$\n",
"is also a solution of the discrete equations.\n",
"[vib:exer:verify:gen:linear](#vib:exer:verify:gen:linear) asks you to carry out\n",
"[Exercise 10: Use linear/quadratic functions for verification](vib_undamped.ipynb#vib:exer:verify:gen:linear) asks you to carry out\n",
"all the details.\n",
"\n",
"### Quadratic solution\n",
"\n",
"Choosing $\\uex = bt^2 + Vt + I$, with $b$ arbitrary,\n",
"Choosing $u_\\text{e} = bt^2 + Vt + I$, with $b$ arbitrary,\n",
"fulfills the initial conditions and\n",
"fits the ODE if $F$ is adjusted properly. The solution also solves\n",
"the discrete equations with linear damping. However, this quadratic\n",
"polynomial in $t$ does not fulfill the discrete equations in case\n",
"of quadratic damping, because the geometric mean used in the approximation\n",
"of this term introduces an error.\n",
"Doing [vib:exer:verify:gen:linear](#vib:exer:verify:gen:linear) will reveal\n",
"Doing [Exercise 10](vib_undamped.ipynb#vib:exer:verify:gen:linear) will reveal\n",
"the details. One can fit $F^n$ in the discrete equations such that\n",
"the quadratic polynomial is reproduced by the numerical method (to\n",
"machine precision).\n",
Expand Down Expand Up @@ -742,7 +754,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"This results in a [moving window following the function](${doc_notes}/vib/html//mov-vib/vib_generalized_dt0.05/index.html) on the screen.\n",
"This results in a [moving window following the function](mov-vib/vib_generalized_dt0.05/index.html) on the screen.\n",
"[Figure](#vib:ode2:fig:demo) shows a part of the time series.\n",
"\n",
"<!-- dom:FIGURE: [fig-vib/vib_gen_demo.png, width=600 frac=1.0] Damped oscillator excited by a sinusoidal function. <div id=\"vib:ode2:fig:demo\"></div> -->\n",
Expand All @@ -758,7 +770,7 @@
"## The Euler-Cromer scheme for the generalized model\n",
"<div id=\"vib:ode2:EulerCromer\"></div>\n",
"\n",
"The ideas of the Euler-Cromer method from the section [vib:model2x2:EulerCromer](#vib:model2x2:EulerCromer)\n",
"The ideas of the [Euler-Cromer method](vib_undamped.ipynb#vib:model2x2:EulerCromer)\n",
"carry over to the generalized model. We write ([1](#vib:ode2))\n",
"as two equations for $u$ and $v=u^{\\prime}$. The first equation is taken as the\n",
"one with $v^{\\prime}$ on the left-hand side:"
Expand Down Expand Up @@ -926,7 +938,8 @@
"## The Stoermer-Verlet algorithm for the generalized model\n",
"<div id=\"vib:model2x2:gen:StormerVerlet\"></div>\n",
"\n",
"We can easily apply the ideas from the section [vib:model2x2:StormerVerlet](#vib:model2x2:StormerVerlet) to\n",
"We can easily apply the ideas from [The Stoermer-Verlet algorithm\n",
"](vib_undamped.ipynb#vib:model2x2:StormerVerlet) section to\n",
"extend that method to the generalized model"
]
},
Expand Down Expand Up @@ -1012,7 +1025,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"mathcal{I}_t is natural to introduce a staggered mesh (see the section [vib:model2x2:staggered](#vib:model2x2:staggered)) and seek $u$ at mesh points $t_n$ (the numerical value is\n",
"It is natural to introduce a staggered mesh (see section [The Euler-Cromer scheme on a staggered mesh](vib_undapmed.ipynb#vib:model2x2:staggered)) and seek $u$ at mesh points $t_n$ (the numerical value is\n",
"denoted by $u^n$) and $v$ between mesh points at $t_{n+1/2}$ (the numerical\n",
"value is denoted by $v^{n+\\frac{1}{2}}$).\n",
"A centered difference approximation to ([26](#vib:ode2:staggered:ueq))-([25](#vib:ode2:staggered:veq)) can then be written in operator notation as"
Expand Down Expand Up @@ -1175,8 +1188,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"The initial conditions are derived at the end of\n",
"the section [vib:model2x2:staggered](#vib:model2x2:staggered):"
"The initial conditions are derived at the end of [The Euler-Cromer scheme on a staggered mesh](vib_undapmed.ipynb#vib:model2x2:staggered) section:"
]
},
{
Expand Down Expand Up @@ -1221,7 +1233,7 @@
"<div id=\"vib:ode2:PEFRL\"></div>\n",
"\n",
"A variant of the Euler-Cromer type of algorithm, which provides an\n",
"error $\\Oof{\\Delta t^4}$ if $f(v)=0$, is called PEFRL\n",
"error $\\mathcal{O}(\\Delta t^4)$ if $f(v)=0$, is called PEFRL\n",
"[[PEFRL_2002]](#PEFRL_2002). This algorithm is very well suited for integrating\n",
"dynamic systems (especially those without damping) over very long time\n",
"periods. Define"
Expand Down Expand Up @@ -1638,7 +1650,7 @@
"source": [
"The advantage of the backward difference is that the damping term is\n",
"evaluated using known values $u^n$ and $u^{n-1}$ only.\n",
"Extend the [`vib.py`](${src_vib}/vib.py) code with a scheme based\n",
"Extend the [`vib.py`](https://github.com/devitocodes/devito_book/blob/master/fdm-devito-notebooks/01_vib/src-vib/vib.py) code with a scheme based\n",
"on using backward differences in the damping terms. Add statements\n",
"to compare the original approach with centered difference and the\n",
"new idea launched in this exercise. Perform numerical experiments\n",
Expand Down Expand Up @@ -2160,13 +2172,8 @@
"$u$ is marched forward by a one-sided Forward Euler scheme applied\n",
"to the first equation, and\n",
"thereafter $v$ can be marched forward by a Backward Euler scheme in the\n",
"second\n",
"% if BOOK == \"book\":\n",
"equation, see in the section [vib:model2x2:EulerCromer](#vib:model2x2:EulerCromer).\n",
"% else:\n",
"equation.\n",
"% endif\n",
"Express the idea in operator notation and write out the\n",
"second <!-- % if BOOK == \"book\": -->\n",
"equation, see in [The Euler-Cromer method](vib_undamped.ipynb#vib:model2x2:EulerCromer) section. <!-- % else: --><!-- equation. --><!-- % endif --> Express the idea in operator notation and write out the\n",
"scheme. Unfortunately, the backward difference for the $v$ equation\n",
"creates a nonlinearity $|v^{n+1}|v^{n+1}$. To linearize this\n",
"nonlinearity, use the known value $v^n$ inside the absolute value\n",
Expand All @@ -2177,16 +2184,19 @@
"and the linearization trick play together in \"the right way\" such that\n",
"the scheme is as good as when we (in the section [A staggered Euler-Cromer scheme for a generalized model](#vib:ode2:staggered))\n",
"carefully apply centered differences and a geometric mean on a\n",
"staggered mesh to achieve second-order accuracy.\n",
"% if BOOK == \"book\":\n",
"There is a\n",
"staggered mesh to achieve second-order accuracy.<!-- % if BOOK == \"book\": --> There is a\n",
"difference in the handling of the initial conditions, though, as\n",
"explained at the end of the section [vib:model2x2:EulerCromer](#vib:model2x2:EulerCromer).\n",
"% endif\n",
"Filename: `vib_gen_bwdamping`.\n",
"explained at the end of [The Euler-Cromer method](vib_undamped.ipynb#vib:model2x2:EulerCromer) section.<!-- % endif --> Filename: `vib_gen_bwdamping`.\n",
"\n",
"<!-- --- end exercise --- -->"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand All @@ -2205,7 +2215,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
"version": "3.8.2"
}
},
"nbformat": 4,
Expand Down