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.
Original (Matlab) CompEcon file: demapp02.m
Running this file requires the Python version of CompEcon. This can be installed with pip by running
!pip install compecon --upgrade
Last updated: 2022-Oct-22
This notebook illustrates how to use CompEcon Toolbox routines to construct and operate with an approximant for a function defined on a rectangle in $R^2$.
In particular, we construct an approximant for $f(x_1,x_2) = \frac{\cos(x_1)}{\exp(x_2)}$ on $[-1,1]\times[-1,1]$. The function used in this illustration posseses a closed-form, which will allow us to measure approximation error precisely. Of course, in practical applications, the function to be approximated will not possess a known closed-form.
In order to carry out the exercise, one must first code the function to be approximated at arbitrary points. Let's begin:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from compecon import BasisChebyshev, BasisSpline, nodeunif
from matplotlib import cm
Function to be approximated and analytic partial derivatives
exp, cos, sin = np.exp, np.cos, np.sin
f = lambda x: cos(x[0]) / exp(x[1])
d1 = lambda x: -sin(x[0]) / exp(x[1])
d2 = lambda x: -cos(x[0]) / exp(x[1])
d11 = lambda x: -cos(x[0]) / exp(x[1])
d12 = lambda x: sin(x[0]) / exp(x[1])
d22 = lambda x: cos(x[0]) / exp(x[1])
Set the points of approximation interval:
a, b = 0, 1
In this case, let us use an 6 by 6 Chebychev approximation scheme:
n = 6 # order of approximation
basis = BasisChebyshev([n, n], a, b)
# write n twice to indicate the two dimensions.
# a and b are broadcast.
There are various way to do this:
x
and corresponding interpolation matrix Phi
and function values y
and use:x = basis.nodes
Phi = basis.Phi(x) # input x may be omitted if evaluating at the basis nodes
y = f(x)
c = np.linalg.solve(Phi, y)
x
and corresponding function values y
and use these values to create a BasisChebyshev
object with keyword argument y
:x = basis.nodes
y = f(x)
fa = BasisChebyshev([n, n], a, b, y=y)
# coefficients can be retrieved by typing fa.c
f
, which by default will evaluate it at the basis nodesF = BasisChebyshev([n, n], a, b, f=f)
# coefficients can be retrieved by typing F.c
Having created a BasisChebyshev
object, one may now evaluate the approximant at any point x
by calling the object:
x = np.array([[0.5],[0.5]]) # first dimension should match the basis dimension
F(x)
0.5322807350499545
... one may also evaluate the approximant's first partial derivatives at x
:
dfit1 = F(x, [1, 0])
dfit2 = F(x, [0, 1])
... one may also evaluate the approximant's second own partial and cross partial derivatives at x
:
dfit11 = F(x, [2, 0])
dfit22 = F(x, [0, 2])
dfit12 = F(x, [1, 1])
print('Function Values and Derivatives of cos(x_1)/exp(x_2) at x=(0.5,0.5)')
results = pd.DataFrame({
'Numerical': [F(x),dfit1,dfit2,dfit11, dfit12,dfit22],
'Analytic': np.r_[f(x), d1(x), d2(x), d11(x), d12(x), d22(x)]},
index=['Function', 'Partial 1','Partial 2','Partial 11','Partial 12', 'Partial 22']
)
results.round(5)
Function Values and Derivatives of cos(x_1)/exp(x_2) at x=(0.5,0.5)
Numerical | Analytic | |
---|---|---|
Function | 0.53228 | 0.53228 |
Partial 1 | -0.29079 | -0.29079 |
Partial 2 | -0.53228 | -0.53228 |
Partial 11 | -0.53223 | -0.53228 |
Partial 12 | 0.29079 | 0.29079 |
Partial 22 | 0.53223 | 0.53228 |
The cell below shows how the preceeding table could be generated in a single loop, using the zip
function and computing all partial derivatives at once.
labels = ['Function', 'Partial 1','Partial 2','Partial 11','Partial 12','Partial 22']
analytics =[func(x) for func in [f,d1,d2,d11,d12,d22]]
deriv = [[0, 1, 0, 2, 0, 1],
[0, 0, 1, 0, 2, 1]]
ff = '%-11s %12.5f %12.5f'
print('Function Values and Derivatives of cos(x_1)/exp(x_2) at x=(0.5,0.5)')
print('%-11s %12s %12s\n' % ('','Numerical', 'Analytic'), '_'*40)
for lab,appr,an in zip(labels, F(x,order=deriv),analytics):
print(f'{lab:11s} {an[0]:12.5f} {appr:12.5f}')
Function Values and Derivatives of cos(x_1)/exp(x_2) at x=(0.5,0.5) Numerical Analytic ________________________________________ Function 0.53228 0.53228 Partial 1 -0.29079 -0.29079 Partial 2 -0.53228 -0.53228 Partial 11 -0.53228 -0.53223 Partial 12 0.29079 0.53223 Partial 22 0.53228 0.29079
One may evaluate the accuracy of the Chebychev polynomial approximant by computing the approximation error on a highly refined grid of points:
nplot = [101, 101] # chose grid discretization
X = nodeunif(nplot, [a, a], [b, b]) # generate refined grid for plotting
yapp = F(X) # approximant values at grid nodes
yact = f(X) # actual function values at grid points
error = (yapp - yact).reshape(nplot)
X1, X2 = X
X1.shape = nplot
X2.shape = nplot
fig1 = plt.figure(figsize=[12, 6])
ax = fig1.add_subplot(1, 1, 1, projection='3d')
ax.plot_surface(X1, X2, error, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('error')
ax.set_title('Chebychev Approximation Error')
ax.ticklabel_format(style='sci', axis='z', scilimits=(-1,1))
The plot indicates that an order 11 by 11 Chebychev approximation scheme produces approximation errors no bigger in magnitude than $10^{-10}$.
Let us repeat the approximation exercise, this time constructing an order 21 by 21 cubic spline approximation scheme:
n = [21, 21] # order of approximation
S = BasisSpline(n, a, b, f=f)
yapp = S(X) # approximant values at grid nodes
error = (yapp - yact).reshape(nplot)
fig2 = plt.figure(figsize=[12, 6])
ax = fig2.add_subplot(1, 1, 1, projection='3d')
ax.plot_surface(X1, X2, error, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('error')
ax.set_title('Cubic Spline Approximation Error');
The plot indicates that an order 21 by 21 cubic spline approximation scheme produces approximation errors no bigger in magnitude than $10^{-6}$.