""" 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()