<jupyter_start><jupyter_text>CompEcon Toolbox:DemQua08Illustrates integration using Simpson's rule<NAME>, PhDThis demo is based on the original Matlab demo accompanying the  Computational Economics and Finance 2001 textbook by <NAME> and <NAME>.Last updated: 2020-Sep-10  Initial tasks<jupyter_code>if 'google.colab' in str(get_ipython()):
    print("This notebook is running on Google Colab. Installing the compecon package.")
    !pip install compecon
from numpy import poly1d,polyfit, linspace, array
from compecon import qnwsimp, demo
import matplotlib.pyplot as plt
n = 1001
xmin, xmax = -1, 1
xwid = xmax-xmin
x = linspace(xmin, xmax, n)
f = poly1d([2.0, -1.0, 0.5, 5.0])
def fitquad(xi):
    newcoef = polyfit(xi, f(xi), 2)
    return poly1d(newcoef)
def plot_simp(n):
    xi, wi = qnwsimp(n+1, xmin, xmax)
    
    fig, ax = plt.subplots()
    ax.plot(x, f(x), linewidth=3)
    
    for k in range(n//2):
        xii = xi[(2*k):(2*k+3)]
        xiii = linspace(xii[0], xii[2], 125)
        p = fitquad(xii)
        ax.fill_between(xiii, p(xiii), alpha=0.35, color='LightSkyBlue')    
        ax.plot(xiii, p(xiii),color='Tab:Orange', linestyle='--')
    
    plt.vlines(xi, 0, f(xi),'k', linestyle=':')
    plt.hlines(0,xmin-0.1, xmax+0.1,'k',linewidth=2)
    plt.xlim(xmin-0.1, xmax+0.1)
    xtl = ['$x_{%d}$' % i for i in range(n+1)]
    xtl[0] += '=a'
    xtl[n] += '=b'
    plt.xticks(xi, xtl)
    plt.yticks([0],['0'])
    plt.legend([r'$f(x)$', f'$\\tilde{{f}}_{n+1}(x)$'])
    return fig
def plot_simp(n):
    xi, wi = qnwsimp(n+1, xmin, xmax)
    
    fig, ax = plt.subplots()
    ax.plot(x, f(x), linewidth=3)
    
    for k in range(n//2):
        xii = xi[(2*k):(2*k+3)]
        xiii = linspace(xii[0], xii[2], 125)
        p = fitquad(xii)
        ax.fill_between(xiii, p(xiii), alpha=0.35, color='LightSkyBlue')    
        ax.plot(xiii, p(xiii), color='Tab:Orange', linestyle='--')
    
    ax.vlines(xi, 0, f(xi), color='Tab:Orange', linestyle=':')
    ax.axhline(0,color='k',linewidth=2)
    
    xtl = [f'$x_{i}$' for i in range(n+1)]
    xtl[0] += '=a'
    xtl[n] += '=b'
    
    ax.set(xlim=[xmin-0.1, xmax+0.1], xticks=xi, xticklabels=xtl,
           yticks=[0], yticklabels=['0'])
    
    plt.legend([r'$f(x)$', f'$\\tilde{{f}}_{n+1}(x)$'])
    return fig

figs = [plot_simp(n) for n in [2, 4, 8]]
#demo.savefig(figs,name='demqua08')<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup