[notebooks] add parametric delay example

notebooks/parametric_delay.ipynb

@@ -0,0 +1,211 @@
+{
+ "cells": [
+  {
+   "cell_type": "raw",
+   "metadata": {},
+   "source": [
+    "This file is part of the pyMOR project (http://www.pymor.org).\n",
+    "Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.\n",
+    "License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import scipy.linalg as spla\n",
+    "import matplotlib as mpl\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from pymor.basic import *"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def H(s, mu):\n",
+    "    tau = mu['tau']\n",
+    "    return np.array([[np.exp(-s) / (tau * s + 1)]])\n",
+    "\n",
+    "def dH(s, mu):\n",
+    "    tau = mu['tau']\n",
+    "    return np.array([[-(tau * s + tau + 1) * np.exp(-s) / (tau * s + 1) ** 2]])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f = ProjectionParameterFunctional('tau', ())\n",
+    "parameter_space = CubicParameterSpace(f.parameter_type, 0.01, 1)\n",
+    "\n",
+    "fom = TransferFunction(NumpyVectorSpace(1), NumpyVectorSpace(1),\n",
+    "                       H, dH,\n",
+    "                       parameter_space=parameter_space)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Bode plot"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mu_list_short = [0.01, 0.1, 1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "w = np.logspace(-2, 4, 100)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "for mu in mu_list_short:\n",
+    "    fom.mag_plot(w, ax=ax, mu=mu, label=fr'$\\tau = {mu}$')\n",
+    "ax.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "w_list = np.logspace(-2, 4, 100)\n",
+    "mu_list = np.logspace(-2, 0, 50)\n",
+    "\n",
+    "fom_w_mu = np.zeros((len(w_list), len(mu_list)))\n",
+    "for i, mu in enumerate(mu_list):\n",
+    "    fom_w_mu[:, i] = spla.norm(fom.bode(w_list, mu=mu), axis=(1, 2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "out = ax.contourf(w_list, mu_list, fom_w_mu.T,\n",
+    "                  norm=mpl.colors.LogNorm(),\n",
+    "                  levels=np.logspace(np.log10(fom_w_mu.min()), np.log10(fom_w_mu.max()), 100))\n",
+    "ax.set_xlabel(r'Frequency $\\omega$')\n",
+    "ax.set_ylabel(r'Parameter $\\mu$')\n",
+    "ax.set_xscale('log')\n",
+    "ax.set_yscale('log')\n",
+    "fig.colorbar(out, ticks=np.logspace(-4, 1, 6))\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# TF-IRKA"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r = 10\n",
+    "roms_tf_irka = []\n",
+    "for mu in mu_list_short:\n",
+    "    tf_irka = TFIRKAReductor(fom, mu=mu)\n",
+    "    rom = tf_irka.reduce(r, conv_crit='h2', maxit=1000, num_prev=5)\n",
+    "    roms_tf_irka.append(rom)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "for mu, rom in zip(mu_list_short, roms_tf_irka):\n",
+    "    poles = rom.poles()\n",
+    "    ax.plot(poles.real, poles.imag, '.', label=fr'$\\tau = {mu}$')\n",
+    "ax.set_title(\"Poles of TF-IRKA's ROMs\")\n",
+    "ax.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "for mu, rom in zip(mu_list_short, roms_tf_irka):\n",
+    "    rom.mag_plot(w, ax=ax, label=fr'$\\tau = {mu}$')\n",
+    "ax.set_title(\"Bode plot of TF-IRKA's ROMs\")\n",
+    "ax.legend()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "for mu, rom in zip(mu_list_short, roms_tf_irka):\n",
+    "    (fom - rom).mag_plot(w, ax=ax, mu=mu, label=fr'$\\tau = {mu}$')\n",
+    "ax.set_title(\"Bode plot of error systems\")\n",
+    "ax.legend()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}