CompEcon Toolbox:
DemApp04
Approximating Runge's function
Randall Romero Aguilar, PhD

This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.

Last updated: 2020-Sep-08

About

Uniform-node and Chebyshev-node polynomial approximation of Runge's function and compute condition numbers of associated interpolation matrices

Initial tasks

In [ ]:
if 'google.colab' in str(get_ipython()):
    print("This notebook is running on Google Colab. Installing the compecon package.")
    !pip install compecon
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm, cond
from compecon import BasisChebyshev, demo
import warnings

warnings.simplefilter('ignore')

Runge function

In [ ]:
runge = lambda x: 1 / (1 + 25 * x ** 2)

Set points of approximation interval

In [ ]:
a, b = -1, 1

Construct plotting grid

In [ ]:
nplot = 1001
x = np.linspace(a, b, nplot)
y = runge(x)

Plot Runge's Function

Initialize data matrices

In [ ]:
n = np.arange(3, 33, 2)
nn = n.size
errunif, errcheb = (np.zeros([nn, nplot]) for k in range(2))
nrmunif, nrmcheb, conunif, concheb = (np.zeros(nn) for k in range(4))

Compute approximation errors on refined grid and interpolation matrix condition numbers

In [ ]:
for i in range(nn):
    # Uniform-node monomial-basis approximant
    xnodes = np.linspace(a, b, n[i])
    c = np.polyfit(xnodes, runge(xnodes), n[i])
    yfit = np.polyval(c, x)
    phi = xnodes.reshape(-1, 1) ** np.arange(n[i])

    errunif[i] = yfit - y
    nrmunif[i] = np.log10(norm(yfit - y, np.inf))
    conunif[i] = np.log10(cond(phi, 2))

    # Chebychev-node Chebychev-basis approximant
    yapprox = BasisChebyshev(n[i], a, b, f=runge)
    yfit = yapprox(x)  # [0] no longer needed?  # index zero is to eliminate one dimension
    phi = yapprox.Phi()
    errcheb[i] = yfit - y
    nrmcheb[i] = np.log10(norm(yfit - y, np.inf))
    concheb[i] = np.log10(cond(phi, 2))

Plot Chebychev- and uniform node polynomial approximation errors

In [ ]:
figs = []
fig1, ax = plt.subplots()
ax.plot(x, y)
ax.text(-0.8, 0.8, r'$y = \frac{1}{1+25x^2}$', fontsize=18)
ax.set(xticks=[], title="Runge's Function", xlabel='', ylabel='y');
figs.append(fig1)
In [ ]:
fig2, ax = plt.subplots()
ax.hlines(0, a, b, 'gray', '--')
ax.plot(x, errcheb[4], label='Chebychev Nodes')
ax.plot(x, errunif[4], label='Uniform Nodes')
ax.legend(loc='upper center')
ax.set(title="Runge's Function $11^{th}$-Degree\nPolynomial Approximation Error.", xlabel='x', ylabel='Error')
figs.append(fig2)

Plot approximation error per degree of approximation

In [ ]:
fig3, ax = plt.subplots()
ax.plot(n, nrmcheb, label='Chebychev Nodes')
ax.plot(n, nrmunif, label='Uniform Nodes')
ax.legend(loc='upper left')
ax.set(title="Log10 Polynomial Approximation Error for Runge's Function",xlabel='', ylabel='Log10 Error', xticks=[])
figs.append(fig3)
In [ ]:
fig4, ax = plt.subplots()
ax.plot(n, concheb, label='Chebychev Polynomial Basis')
ax.plot(n, conunif, label='Mononomial Basis')
ax.legend(loc='upper left')
ax.set(title="Log10 Interpolation Matrix Condition Number",
       xlabel='Degree of Approximating Polynomial',
       ylabel='Log10 Condition Number')
figs.append(fig4)

Save all figures to disc

In [ ]:
#demo.savefig(figs, name='demapp04')