"""
fit_cubic.py
OVERVIEW
This script demonstrates how to solve a small system of nonlinear equations
numerically using Python and SciPy, as well as how to perform basic operations
on polynomials. The script does the following:
(1) It uses one of SciPy's root solvers to find the cubic polynomial p(x) such
that the polynomial and its derivative q(x) have specified values for two values
of x.
(2) It plots the resulting polynomial.
To modify the inputs, edit Section 2 below.
AUTHOR
Phillip M. Feldman
"""
# Section 1: Import from modules.
import pdb
from matplotlib import pyplot
from numpy import linspace, polyder, polyval
from scipy import optimize
# Section 2: Define x_0, x_1, and the corresponding values (p's) and derivatives
# (q's) at those two locations:
x_0 = 0.0
x_1 = 1.0
p_x_0= 0.0
q_x_0= 1.0
p_x_1= 0.0
q_x_1= 1.0
# Section 3: Define error function `E`. This function returns a length-4 list
# of errors (differences) between each function value and derivative and the
# corresponding target value. Like `numpy.polyval`, the function expects
# coefficients to be ordered starting with the highest power.
def E(polynomial):
# Calculate derivative polynomial:
derivative= polyder(polynomial)
return [
polyval(polynomial, x_0) - p_x_0,
polyval(derivative, x_0) - q_x_0,
polyval(polynomial, x_1) - p_x_1,
polyval(derivative, x_1) - q_x_1,
]
# Section 4: Run solver to get polynomial coefficients.
coefs= optimize.broyden1(E, [0, 0, 0, 1], f_tol=1e-14)
# Section 5: Plot the results.
# Create figure window:
fig= pyplot.figure(figsize=[9,7], dpi=120, facecolor=[1,1,1])
x= linspace(0.0, 1.0, 50)
pyplot.plot(x, polyval(coefs, x), linewidth=2)
pyplot.grid(True)
pyplot.show()
# pdb.set_trace()