Compute root of f(x)=\exp(-x)-1
Contents
Compute root of \(f(x)=\exp(-x)-1\)¶
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: demslv01.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-Sept-04
About¶
Compute root of \(f(x)=\exp(-x)-1\) using Newton and secant methods. Initial value generated randomly. True root is \(x=0\).
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from compecon import NLP, tic, toc
Set up the problem¶
def f(x):
fval = np.exp(-x) - 1
fjac = -np.exp(-x)
return fval, fjac
problem = NLP(f, all_x=True)
Randomly generate starting point¶
np.random.RandomState(seed=2022)
problem.x0 = -10 * np.random.rand(1)
Compute root using Newton method¶
t0 = tic()
x1 = problem.newton()
t1 = 100 * toc(t0)
n1, x_newton = problem.fnorm, problem.x_sequence
x_newton
| x_0 | |
|---|---|
| iteration | |
| 0 | -9.364871e+00 |
| 1 | -8.364957e+00 |
| 2 | -7.365189e+00 |
| 3 | -6.365822e+00 |
| 4 | -5.367542e+00 |
| 5 | -4.372207e+00 |
| 6 | -3.384831e+00 |
| 7 | -2.418714e+00 |
| 8 | -1.507750e+00 |
| 9 | -7.291576e-01 |
| 10 | -2.114727e-01 |
| 11 | -2.086409e-02 |
| 12 | -2.161492e-04 |
| 13 | -2.335855e-08 |
| 14 | -3.674602e-16 |
Compute root using Broyden method¶
t0 = tic()
x2 = problem.broyden()
t2 = 100 * toc(t0)
n2, x_broyden = problem.fnorm, problem.x_sequence
x_broyden
| x_0 | |
|---|---|
| iteration | |
| 0 | -9.364871e+00 |
| 1 | -8.364957e+00 |
| 2 | -7.783086e+00 |
| 3 | -7.046272e+00 |
| 4 | -6.370431e+00 |
| 5 | -5.671773e+00 |
| 6 | -4.983130e+00 |
| 7 | -4.293003e+00 |
| 8 | -3.608176e+00 |
| 9 | -2.930681e+00 |
| 10 | -2.268779e+00 |
| 11 | -1.636439e+00 |
| 12 | -1.059074e+00 |
| 13 | -5.763789e-01 |
| 14 | -2.355628e-01 |
| 15 | -5.942547e-02 |
| 16 | -6.663203e-03 |
| 17 | -1.958079e-04 |
| 18 | -6.516081e-07 |
| 19 | -6.379284e-11 |
Print results¶
print('Hundredths of seconds required to compute root of exp(-x)-1,')
print('via Newton and Broyden methods, starting at x = %4.2f.' % problem.x0)
pd.DataFrame({
'Time': [t1, t2],
'Norm of f': [n1, n2],
'Final x': [x1, x2]},
index=['Newton', 'Broyden']
)
Hundredths of seconds required to compute root of exp(-x)-1,
via Newton and Broyden methods, starting at x = -9.36.
| Time | Norm of f | Final x | |
|---|---|---|---|
| Newton | 0.299716 | 4.440892e-16 | -3.674602e-16 |
| Broyden | 0.747013 | 6.379275e-11 | -6.379284e-11 |
View current options for solver¶
print(problem.opts)
Options for solving a NLP:
all_x = True
initb = None
initi = False
maxit = 100
maxsteps = 10
method = broyden
show = False
squeeze = True
tol = 1.4901161193847656e-08
transform = ssmooth
Describe the options¶
print(problem.opts.__doc__)
A container for options to find a zero for a NLP or MCP
Attributes: default in brackets
method: either ['newton'] or 'broyden'
maxit: maximum number of iterations [100]
maxsteps: maximum number of backsteps [10]
tol: convergence tolerance [sqrt(machine eps)]
initb: an initial inverse Jacobian aprroximation matrix [None]
initi: if True, use the identity matrix to initialize Jacobian,
if [False], a numerical Jacobian will be used
transform: either ['ssmooth'] or 'minmax', required for MC problems
show: print to screen if [True], quiet if False
all_x: whether to output the full solution sequence too [False]
squeeze: if problem has ony one dimension, return solution as scalar [True]
Plot the convergence¶
b = -abs(problem.x0)
a = -b
xx = np.linspace(a, b, 100)
fig, ax = plt.subplots()
ax.hlines(0, a, b, 'gray')
ax.plot(xx, f(xx)[0], 'b-', alpha=0.4)
ax.plot(x_newton,f(x_newton)[0],'ro:', alpha=0.4);
=exp(-x)-1_18_0.png)