{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One key concept in physics is the conservation of energy---in an isolated system, the total energy cannot increase or decrease, it simply remains constant. Thus, the simple lossless harmonic oscillator model defined by the equation\n",
"\n",
"$$ M\\frac{d^2 u}{dt^2} = -K u $$\n",
"\n",
"has an associated energy $H$ defined as\n",
"\n",
"$$ H(t) = \\frac{M}{2} \\left( \\frac{du}{dt} \\right)^2 + \\frac{K}{2} u^2, \\qquad \\frac{d}{dt} H =0,$$\n",
"\n",
"which is conserved. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Such a concept can also be applied in a numerical simulation, where the existence of a numerical energy-like quantity with certain properties can become a powerful tool in guaranteeing the stability of the scheme. \n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Discrete operators\n",
"\n",
"A finite difference scheme for the lumped oscillator equation above has been previously discussed [here](http://blog.albertotorin.it/blog/anatomy-fd-sho/). Before jumping into the analysis of numerical energy, it is worth introducing some useful notation.\n",
"\n",
"Given a discrete variable $u^n$ defined at integer multiples $n$ of a time step $k$, the following discrete operators can be defined:\n",
"\n",
"$$ \\delta_{tt} u^n = \\frac{u^{n+1} - 2u^n + u^{n-1}}{k^2}, \\qquad \\delta_{t+} u^n = \\frac{u^{n+1} - u^n}{k}, \\qquad \\delta_{t-} u^n = \\frac{u^{n} - u^{n-1}}{k}. $$\n",
"\n",
"These are just a shorthand for the second and first order difference operators. With this notation, the finite difference scheme for the SHO can be written in a much more compact way as\n",
"\n",
"$$ M\\delta_{tt} u^n = - K u^n. $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Discrete energy\n",
"\n",
"It is possible to show that the above scheme as an associated discrete energy quantity defined as\n",
"\n",
"$$ H^{n+1/2} = \\frac{M}{2} \\left( \\delta_{t+} u^n \\right)^2 + \\frac{K}{2} u^n \\, u^{n+1} $$\n",
"\n",
"which is conserved, such that\n",
"\n",
"$$ \\delta_{t-} H^{n+1/2} = 0. $$\n",
"\n",
"Notice that the discrete energy is defined at interleaved steps to the displacement $u^n$. The conservation of energy in the discrete case holds assuming infinite precision arithmetic, but breaks down when finite precision representation of numbers (such as floating point) is used. Variations, however, are small, and confined to the last digits of the values, as will be seen shortly."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numerical simulation\n",
"\n",
"The calculation of the discrete energy can be easily implemented in the code given [here](https://github.com/atorin/FiniteDifferenceSchemes/blob/master/lossy_sho.jl). All it is necessary to do is to define a function that calculates the new value of $H^{n+1/2}$ at every time step. In order to simplify the implementation, it is possible to scale the value by $M$ and introduce the parameter $\\omega_0^2 = K/M.$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"H_calc (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function H_calc(u, u1, k, ω0)\n",
" v = (u-u1)/k\n",
" H = 0.5 * v^2 + 0.5 * ω0^2 * u * u1\n",
" return H\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The finite difference code for the SHO with energy is the following."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"sho (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function sho(f)\n",
"\n",
"# Initial conditions\n",
"x0=0.01\n",
"v0=0\n",
"\n",
"SR=44100 # Sampling rate\n",
"T=1 # Duration of the simulation (s)\n",
"\n",
"# Calculate the time step\n",
"k=1/SR\n",
"# Calculate the number of iterations\n",
"NF=floor(SR*T)\n",
"\n",
"# Calculate ω0. This will be useful later\n",
"ω0=2*pi*f\n",
"\n",
"# Check the stability condition\n",
"if k>(2/ω0)\n",
" println(\"Stability condition violated!\")\n",
" return\n",
"end\n",
"\n",
"# Initialise and allocate displacement\n",
"u0=x0\n",
"u1=x0+k*v0\n",
"out=zeros(NF)\n",
"out[1]=u0\n",
"out[2]=u1\n",
" \n",
"# Allocate energy vector\n",
"H = zeros(NF-1)\n",
"H[1] = H_calc(u1, u0, k, ω0)\n",
"\n",
"# Start the Main Loop\n",
"for nn=3:NF\n",
" # Update the recursion\n",
" u = 2*u1 - ω0^2*k^2*u1 - u0\n",
"\n",
" # Store the new value for the position\n",
" out[nn]=u\n",
" \n",
" # Store the new value for the energy\n",
" H[nn-1] = H_calc(u, u1, k, ω0)\n",
" \n",
" # Reset the known values of the recursion\n",
" u0, u1 = u1, u;\n",
"\n",
"end\n",
"\n",
" # At the end of the loop, return output\n",
" return out, H\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Energy conservation\n",
"\n",
"It is possible now to run a simulation and check if the discrete energy is conserved... "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"