nbconvert parametric_* notebooks to demos

all other demos already had a demo twin
parent bbeabc63
#!/usr/bin/env python
# coding: utf-8
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
# In[ ]:
import numpy as np
import scipy.linalg as spla
import matplotlib as mpl
import matplotlib.pyplot as plt
from pymor.basic import *
# # Model
# In[ ]:
def H(s, mu):
tau = mu['tau'][0]
return np.array([[np.exp(-s) / (tau * s + 1)]])
def dH(s, mu):
tau = mu['tau'][0]
return np.array([[-(tau * s + tau + 1) * np.exp(-s) / (tau * s + 1) ** 2]])
# In[ ]:
fom = TransferFunction(1, 1,
H, dH,
parameters={'tau': 1})
# # Magnitude plot
# In[ ]:
mu_list_short = [0.01, 0.1, 1]
# In[ ]:
w = np.logspace(-2, 4, 100)
fig, ax = plt.subplots()
for mu in mu_list_short:
fom.mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.legend()
plt.show()
# In[ ]:
w_list = np.logspace(-2, 4, 100)
mu_list = np.logspace(-2, 0, 50)
fom_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
fom_w_mu[:, i] = spla.norm(fom.freq_resp(w_list, mu=mu), axis=(1, 2))
# In[ ]:
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, fom_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(np.log10(fom_w_mu.min()), np.log10(fom_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-4, 1, 6))
plt.show()
# # TF-IRKA
# In[ ]:
r = 10
roms_tf_irka = []
for mu in mu_list_short:
tf_irka = TFIRKAReductor(fom, mu=mu)
rom = tf_irka.reduce(r, conv_crit='h2', maxit=1000, num_prev=5)
roms_tf_irka.append(rom)
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
poles = rom.poles()
ax.plot(poles.real, poles.imag, '.', label=fr'$\tau = {mu}$')
ax.set_title("Poles of TF-IRKA's ROMs")
ax.legend()
plt.show()
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
rom.mag_plot(w, ax=ax, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of TF-IRKA's ROMs")
ax.legend()
plt.show()
# In[ ]:
fig, ax = plt.subplots()
for mu, rom in zip(mu_list_short, roms_tf_irka):
(fom - rom).mag_plot(w, ax=ax, mu=mu, label=fr'$\tau = {mu}$')
ax.set_title("Magnitude plot of error systems")
ax.legend()
plt.show()
#!/usr/bin/env python
# coding: utf-8
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
# In[ ]:
# In[ ]:
import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt
import matplotlib as mpl
from pymor.basic import *
from pymor.core.config import config
from pymor.core.logger import set_log_levels
set_log_levels({'pymor.algorithms.gram_schmidt.gram_schmidt': 'WARNING'})
set_defaults({'pymor.discretizers.builtin.gui.jupyter.get_visualizer.backend': 'not pythreejs'})
# # Model
# In[ ]:
p = InstationaryProblem(
StationaryProblem(
domain=LineDomain([0.,1.], left='robin', right='robin'),
diffusion=LincombFunction([ExpressionFunction('(x[...,0] <= 0.5) * 1.', 1),
ExpressionFunction('(0.5 < x[...,0]) * 1.', 1)],
[1,
ProjectionParameterFunctional('diffusion')]),
robin_data=(ConstantFunction(1., 1), ExpressionFunction('(x[...,0] < 1e-10) * 1.', 1)),
outputs=(('l2_boundary', ExpressionFunction('(x[...,0] > (1 - 1e-10)) * 1.', 1)),),
),
ConstantFunction(0., 1),
T=3.
)
fom, _ = discretize_instationary_cg(p, diameter=1/100, nt=100)
# In[ ]:
fom.visualize(fom.solve(mu=0.1))
# In[ ]:
fom.visualize(fom.solve(mu=1))
# In[ ]:
fom.visualize(fom.solve(mu=10))
# In[ ]:
lti = fom.to_lti()
# # System analysis
# In[ ]:
print(f'order of the model = {lti.order}')
print(f'number of inputs = {lti.dim_input}')
print(f'number of outputs = {lti.dim_output}')
# In[ ]:
mu_list = [0.1, 1, 10]
fig, ax = plt.subplots(len(mu_list), 1, sharex=True, sharey=True)
for i, mu in enumerate(mu_list):
poles = lti.poles(mu=mu)
ax[i].plot(poles.real, poles.imag, '.')
ax[i].set_xscale('symlog')
ax[i].set_title(fr'$\mu = {mu}$')
fig.suptitle('System poles')
fig.subplots_adjust(hspace=0.5)
plt.show()
# In[ ]:
mu_list = [0.1, 1, 10]
fig, ax = plt.subplots()
w = np.logspace(-1, 3, 100)
for mu in mu_list:
lti.mag_plot(w, ax=ax, mu=mu, label=fr'$\mu = {mu}$')
ax.legend()
plt.show()
# In[ ]:
w_list = np.logspace(-1, 3, 100)
mu_list = np.logspace(-1, 1, 20)
lti_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
lti_w_mu[:, i] = spla.norm(lti.freq_resp(w, mu=mu), axis=(1, 2))
# In[ ]:
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, lti_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(-16, np.log10(lti_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-16, 0, 17))
plt.show()
# In[ ]:
mu_list = [0.1, 1, 10]
fig, ax = plt.subplots()
for mu in mu_list:
hsv = lti.hsv(mu=mu)
ax.semilogy(range(1, len(hsv) + 1), hsv, label=fr'$\mu = {mu}$')
ax.set_title('Hankel singular values')
ax.legend()
plt.show()
# In[ ]:
fig, ax = plt.subplots()
mu_fine = np.logspace(-1, 1, 20)
h2_norm_mu = [lti.h2_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, h2_norm_mu, label=r'$\mathcal{H}_2$-norm')
if config.HAVE_SLYCOT:
hinf_norm_mu = [lti.hinf_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hinf_norm_mu, label=r'$\mathcal{H}_\infty$-norm')
hankel_norm_mu = [lti.hankel_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hankel_norm_mu, label='Hankel norm')
ax.set_xlabel(r'$\mu$')
ax.set_title('System norms')
ax.legend()
plt.show()
# # Balanced truncation
# In[ ]:
def reduction_errors(lti, r, mu_fine, reductor, **kwargs):
h2_err_mu = []
hinf_err_mu = []
hankel_err_mu = []
for mu in mu_fine:
rom_mu = reductor(lti, mu=mu, **kwargs).reduce(r)
h2_err_mu.append((lti - rom_mu).h2_norm(mu=mu) / lti.h2_norm(mu=mu))
if config.HAVE_SLYCOT:
hinf_err_mu.append((lti - rom_mu).hinf_norm(mu=mu) / lti.hinf_norm(mu=mu))
hankel_err_mu.append((lti - rom_mu).hankel_norm(mu=mu) / lti.hankel_norm(mu=mu))
return h2_err_mu, hinf_err_mu, hankel_err_mu
# In[ ]:
r = 5
mu_fine = np.logspace(-1, 1, 10)
h2_bt_err_mu, hinf_bt_err_mu, hankel_bt_err_mu = reduction_errors(lti, r, mu_fine, BTReductor)
# In[ ]:
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_bt_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_bt_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_bt_err_mu, '.-', label='Hankel')
ax.set_xlabel(r'$\mu$')
ax.set_title('Balanced truncation errors')
ax.legend()
plt.show()
# # Iterative Rational Krylov Algorithm (IRKA)
# In[ ]:
h2_irka_err_mu, hinf_irka_err_mu, hankel_irka_err_mu = reduction_errors(lti, r, mu_fine, IRKAReductor)
# In[ ]:
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_irka_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_irka_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_irka_err_mu, '.-', label='Hankel')
ax.set_xlabel(r'$\mu$')
ax.set_title('IRKA errors')
ax.legend()
plt.show()
# # Two-Sided Iteration Algorithm (TSIA)
# In[ ]:
h2_tsia_err_mu, hinf_tsia_err_mu, hankel_tsia_err_mu = reduction_errors(lti, r, mu_fine, TSIAReductor)
# In[ ]:
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_tsia_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_tsia_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_tsia_err_mu, '.-', label='Hankel')
ax.set_xlabel(r'$\mu$')
ax.set_title('TSIA errors')
ax.legend()
plt.show()
# # One-sided IRKA
# In[ ]:
h2_osirka_err_mu, hinf_osirka_err_mu, hankel_osirka_err_mu = reduction_errors(
lti, r, mu_fine, OneSidedIRKAReductor, version='V'
)
# In[ ]:
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_osirka_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_osirka_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_osirka_err_mu, '.-', label='Hankel')
ax.set_xlabel(r'$\mu$')
ax.set_title('One-sided IRKA errors')
ax.legend()
plt.show()
This diff is collapsed.
#!/usr/bin/env python
# coding: utf-8
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
# In[ ]:
import numpy as np
import scipy.linalg as spla
import scipy.sparse as sps
import matplotlib.pyplot as plt
import matplotlib as mpl
from pymor.core.config import config
from pymor.models.iosys import LTIModel
from pymor.operators.numpy import NumpyMatrixOperator
from pymor.parameters.functionals import ProjectionParameterFunctional
from pymor.reductors.bt import BTReductor
from pymor.reductors.h2 import IRKAReductor
# # Model
#
# https://morwiki.mpi-magdeburg.mpg.de/morwiki/index.php/Synthetic_parametric_model
# In[ ]:
n = 100 # order of the resulting system
# set coefficients
a = -np.linspace(1e1, 1e3, n // 2)
b = np.linspace(1e1, 1e3, n // 2)
c = np.ones(n // 2)
d = np.zeros(n // 2)
# build 2x2 submatrices
aa = np.empty(n)
aa[::2] = a
aa[1::2] = a
bb = np.zeros(n)
bb[::2] = b
# set up system matrices
Amu = sps.diags(aa, format='csc')
A0 = sps.diags([bb, -bb], [1, -1], shape=(n, n), format='csc')
B = np.zeros((n, 1))
B[::2, 0] = 2
C = np.empty((1, n))
C[0, ::2] = c
C[0, 1::2] = d
# In[ ]:
A0 = NumpyMatrixOperator(A0)
Amu = NumpyMatrixOperator(Amu)
B = NumpyMatrixOperator(B)
C = NumpyMatrixOperator(C)
# In[ ]:
A = A0 + Amu * ProjectionParameterFunctional('mu')
# In[ ]:
lti = LTIModel(A, B, C)
# # Magnitude plot
# In[ ]:
mu_list_short = [1/50, 1/20, 1/10, 1/5, 1/2, 1]
# In[ ]:
w = np.logspace(0.5, 3.5, 200)
fig, ax = plt.subplots()
for mu in mu_list_short:
lti.mag_plot(w, ax=ax, mu=mu, label=fr'$\mu = {mu}$')
ax.legend()
plt.show()
# In[ ]:
w_list = np.logspace(0.5, 3.5, 200)
mu_list = np.linspace(1/50, 1, 50)
lti_w_mu = np.zeros((len(w_list), len(mu_list)))
for i, mu in enumerate(mu_list):
lti_w_mu[:, i] = spla.norm(lti.freq_resp(w_list, mu=mu), axis=(1, 2))
# In[ ]:
fig, ax = plt.subplots()
out = ax.contourf(w_list, mu_list, lti_w_mu.T,
norm=mpl.colors.LogNorm(),
levels=np.logspace(np.log10(lti_w_mu.min()), np.log10(lti_w_mu.max()), 100))
ax.set_xlabel(r'Frequency $\omega$')
ax.set_ylabel(r'Parameter $\mu$')
ax.set_xscale('log')
#ax.set_yscale('log')
fig.colorbar(out, ticks=np.logspace(-2, 1, 7))
plt.show()
# # Hankel singular values
# In[ ]:
fig, ax = plt.subplots()
for mu in mu_list_short:
hsv = lti.hsv(mu=mu)
ax.semilogy(range(1, len(hsv) + 1), hsv, '.-', label=fr'$\mu = {mu}$')
ax.set_title('Hankel singular values')
ax.legend()
plt.show()
# # System norms
# In[ ]:
fig, ax = plt.subplots()
mu_fine = np.linspace(1/50, 1, 20)
h2_norm_mu = [lti.h2_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, h2_norm_mu, '.-', label=r'$\mathcal{H}_2$-norm')
if config.HAVE_SLYCOT:
hinf_norm_mu = [lti.hinf_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hinf_norm_mu, '.-', label=r'$\mathcal{H}_\infty$-norm')
hankel_norm_mu = [lti.hankel_norm(mu=mu) for mu in mu_fine]
ax.plot(mu_fine, hankel_norm_mu, '.-', label='Hankel norm')
ax.set_xlabel(r'$\mu$')
ax.set_title('System norms')
ax.legend()
plt.show()
# # Balanced truncation
# In[ ]:
def reduction_errors(lti, r, mu_fine, method):
h2_err_mu = []
hinf_err_mu = []
hankel_err_mu = []
for mu in mu_fine:
rom_mu = method(lti, r, mu=mu)
h2_err_mu.append((lti - rom_mu).h2_norm(mu=mu) / lti.h2_norm(mu=mu))
if config.HAVE_SLYCOT:
hinf_err_mu.append((lti - rom_mu).hinf_norm(mu=mu) / lti.hinf_norm(mu=mu))
hankel_err_mu.append((lti - rom_mu).hankel_norm(mu=mu) / lti.hankel_norm(mu=mu))
return h2_err_mu, hinf_err_mu, hankel_err_mu
# In[ ]:
r = 20
mu_fine = np.linspace(1/50, 1, 10)
(
h2_bt_err_mu,
hinf_bt_err_mu,
hankel_bt_err_mu
) = reduction_errors(lti, r, mu_fine,
lambda lti, r, mu=None: BTReductor(lti, mu=mu).reduce(r))
# In[ ]:
fig, ax = plt.subplots()
ax.semilogy(mu_fine, h2_bt_err_mu, '.-', label=r'$\mathcal{H}_2$')
if config.HAVE_SLYCOT:
ax.semilogy(mu_fine, hinf_bt_err_mu, '.-', label=r'$\mathcal{H}_\infty$')
ax.semilogy(mu_fine, hankel_bt_err_mu, '.-', label='Hankel')
ax.set_xlabel(r'$\mu$')
ax.set_title('Balanced truncation errors')
ax.legend()
plt.show()
# # Iterative Rational Krylov Algorithm (IRKA)
# In[ ]:
(
h2_irka_err_mu,
hinf_irka_err_mu,
hankel_irka_err_mu
) = reduction_errors(lti, r, mu_fine,
lambda lti, r, mu=mu: IRKAReductor(lti, mu=mu).reduce(r, conv_crit='h2'))
# In[ ]: