182 lines
No EOL
9.2 KiB
Python
182 lines
No EOL
9.2 KiB
Python
# -*- coding: utf-8 -*-
|
|
import numpy as np
|
|
|
|
def cubic_poly(x, a0, a1, a2, a3):
|
|
return a0 + a1*x + a2*x**2 + a3*x**3
|
|
|
|
def polyfit4(x1, x2, x3, x4, y1, y2, y3, y4):
|
|
''' Function to compute the coefficients of the cubic Lagrange polynomial that interpolates the four points (x1, y1), (x2, y2), (x3, y3), (x4, y4).
|
|
The coefficients are computed using a linear least-squares fit to the four points.
|
|
The condition number of the matrix (A^T * A) is critical to the numerical stability of the solution.
|
|
It happens that the minimum condition number is achieved when the four points are equally spaced and span the interval [-pi/3, pi/3].
|
|
'''
|
|
A = np.array([[1, x1, x1**2, x1**3], [1, x2, x2**2, x2**3], [1, x3, x3**2, x3**3], [1, x4, x4**2, x4**3]])
|
|
y = np.array([y1, y2, y3, y4])
|
|
a = np.linalg.solve(A.transpose()@A, A.transpose()@y)
|
|
return a
|
|
|
|
def polyfit4_mpi3_pi3(y1, y2, y3, y4):
|
|
''' Function to compute the coefficients of the cubic Lagrange polynomial that interpolates the four points (-pi/3, y1), (-pi/9, y2), (pi/9, y3), (pi/3, y4).
|
|
The coefficients are computed using a linear least-squares fit to the four points.
|
|
The matrix (A^T * A) is precomputed for efficiency
|
|
'''
|
|
y = np.array([y1, y2, y3, y4])
|
|
AT = np.array([[1, 1, 1, 1],
|
|
[-np.pi/3, -np.pi/9, np.pi/9, np.pi/3],
|
|
[np.pi**2/9, np.pi**2/81, np.pi**2/81, np.pi**2/9],
|
|
[-np.pi**3/27, -np.pi**3/729, np.pi**3/729, np.pi**3/27]]) # A^T
|
|
ATAinv = np.array([[41/64, 0, -405/(64*np.pi**2), 0],
|
|
[0, 3285/(64*np.pi**2), 0, -29889/(64*np.pi**4)],
|
|
[-405/(64*np.pi**2), 0, 6561/(64*np.pi**4), 0],
|
|
[0, -29889/(64*np.pi**4), 0, 295245/(64*np.pi**6)]]) # (A^T * A)^-1
|
|
return ATAinv@(AT@y)
|
|
|
|
def map_interval(x, a1, b1, a2, b2):
|
|
''' Maps a value or array "x" from the interval [a1, b1] to the interval [a2, b2]. '''
|
|
return (x - a1) * (b2 - a2) / (b1 - a1) + a2
|
|
|
|
def cubic_lagrange_minimize(f, a, b, tol=1e-6, callback=None):
|
|
''' Function to find minimum of f over interval [a, b] using cubic Lagrange polynomial interpolation.
|
|
If the function is monotonic, then the minimum is one of the bounds of the interval, and the minimum is found in two iterations.
|
|
The best estimate of the minimum is returned.
|
|
The number of function evaluations is 2 + 2*Niter.
|
|
callback is a function that is called at each iteration with the current estimate of the minimum and other internal parameters of the algorithm.
|
|
The callback function should have the following signature: callback(x_sol, y_sol, i, quadratic_solution, x, a), where quadratic_solution = (delta, x_sol_1, x_sol_2, y_sol_1, y_sol_2) is the solution of the quadratic equation of the first derivative of the Lagrange polynomial, x = (x0, x1, x2, x3) are the endpoints of the interval, and a = (a0, a1, a2, a3) are the coefficients of the Lagrange polynomial.
|
|
'''
|
|
# initialize interval endpoints and function values
|
|
x0, x1, x2, x3 = a, a*2./3. + b*1./3., a*1./3. + b*2./3., b
|
|
f0, f3 = f(x0), f(x3)
|
|
x_prev = (x0 + x3)/2.
|
|
small_coefficient = 1e-9
|
|
Niter = int(np.ceil(np.log(np.abs(b-a)/tol)/np.log(3)))# reduction factor of interval size at each iteration is 3 because we sample two points in the interval each iteration
|
|
|
|
for i in range(Niter):
|
|
# Compute function values at two points in the interval
|
|
x1, x2 = x0*2./3. + x3*1./3., x0*1./3. + x3*2./3.
|
|
f1, f2 = f(x1), f(x2)
|
|
|
|
# Remap the x values to the interval [-pi/3, pi/3] to minimize the condition number of the matrix (A^T * A) in the least-squares fit
|
|
x0m, x1m, x2m, x3m = map_interval(x0, x0, x3, -np.pi/3, np.pi/3), map_interval(x1, x0, x3, -np.pi/3, np.pi/3), map_interval(x2, x0, x3, -np.pi/3, np.pi/3), map_interval(x3, x0, x3, -np.pi/3, np.pi/3)
|
|
|
|
# compute Lagrange polynomial using least-squares fit to 4 points, which is equivalent to the cubic Lagrange polynomial
|
|
a0, a1, a2, a3 = polyfit4_mpi3_pi3(f0, f1, f2, f3) # equivalent to a0, a1, a2, a3 = polyfit4(x0m, x1m, x2m, x3m, f0, f1, f2, f3)
|
|
x_sol_1, x_sol_2, y_sol_1, y_sol_2, delta = None, None, None, None, None
|
|
|
|
# Solve the first derivative of the Lagrange polynomial for a zero
|
|
if np.abs(a3) > small_coefficient:
|
|
delta = -3*a1*a3 + a2**2
|
|
|
|
if delta < 0:
|
|
x_sol = x0m if f0 < f3 else x3m # just choose the interval tha contains the minimum of the linear polynomial
|
|
y_sol = cubic_poly(x_sol, a0, a1, a2, a3)
|
|
else: # solve for the two solutions of the quadratic equation of the first derivative of the Lagrange polynomial
|
|
x_sol_1 = (-a2 + np.sqrt(delta))/(3*a3)
|
|
x_sol_2 = (-a2 - np.sqrt(delta))/(3*a3)
|
|
|
|
y_sol_1 = cubic_poly(x_sol_1, a0, a1, a2, a3)
|
|
y_sol_2 = cubic_poly(x_sol_2, a0, a1, a2, a3)
|
|
|
|
x_sol = x_sol_1 if y_sol_1 < y_sol_2 else x_sol_2
|
|
y_sol = y_sol_1 if y_sol_1 < y_sol_2 else y_sol_2
|
|
|
|
elif np.abs(a2) > small_coefficient: # if a3 is zero, then the Lagrange polynomial is a quadratic polynomial
|
|
x_sol = -a1/(2*a2)
|
|
y_sol = cubic_poly(x_sol, a0, a1, a2, a3)
|
|
else: # if a3 and a2 are zero, then the Lagrange polynomial is a linear polynomial
|
|
x_sol = x0m if f0 < f3 else x3m # just choose the interval tha contains the minimum of the linear polynomial
|
|
y_sol = cubic_poly(x_sol, a0, a1, a2, a3)
|
|
|
|
# transform the solution back to the original interval
|
|
x_sol = map_interval(x_sol, -np.pi/3, np.pi/3, x0, x3)
|
|
if x_sol_1 is not None:
|
|
x_sol_1 = map_interval(x_sol_1, -np.pi/3, np.pi/3, x0, x3)
|
|
if x_sol_2 is not None:
|
|
x_sol_2 = map_interval(x_sol_2, -np.pi/3, np.pi/3, x0, x3)
|
|
|
|
if callback is not None:
|
|
callback(x_sol, y_sol, i, (delta, x_sol_1, x_sol_2, y_sol_1, y_sol_2), (x0, x1, x2, x3), (a0, a1, a2, a3))
|
|
|
|
if np.abs(x_sol - x_prev) < tol:
|
|
break
|
|
|
|
# Determine which interval contains the minimum of the cubic polynomial
|
|
if x_sol < x1:
|
|
x3, f3 = x1, f1
|
|
elif x_sol < x2:
|
|
x0, f0 = x1, f1
|
|
x3, f3 = x2, f2
|
|
else:
|
|
x0, f0 = x2, f2
|
|
|
|
x_prev = x_sol
|
|
|
|
# return best estimate of minimum
|
|
if y_sol < f0 and y_sol < f3:
|
|
return x_sol
|
|
elif f0 < f3:
|
|
return x0
|
|
else:
|
|
return x3
|
|
|
|
def cubic_lagrange_minimize_callback_simple(x_sol, y_sol, i, quadratic_solution, x, a):
|
|
print(f'Iteration {i} : f({x_sol}) = {y_sol}')
|
|
|
|
def cubic_lagrange_minimize_callback_detailed(x_sol, y_sol, i, quadratic_solution, x, a):
|
|
print('--------------------------------------------------------------------')
|
|
print(f'Iteration {i}')
|
|
print('X values : ', x)
|
|
print('Cubic poly coeffs : ', a)
|
|
print(f'Quadratic solution : delta = {quadratic_solution[0]}, x_sol_1 = {quadratic_solution[1]}, x_sol_2 = {quadratic_solution[2]}, y_sol_1 = {quadratic_solution[3]}, y_sol_2 = {quadratic_solution[4]}')
|
|
print(f'Current solution : f({x_sol}) = {y_sol}')
|
|
|
|
if __name__ == '__main__':
|
|
import colorsys
|
|
import matplotlib.pyplot as plt
|
|
|
|
functions_list = [
|
|
lambda x: x**2 + np.sin(5*x),
|
|
lambda x: (x-1)**2,
|
|
lambda x: -x,
|
|
lambda x: np.exp(x),
|
|
lambda x: np.exp(-x),
|
|
lambda x: -np.exp(-x**2),
|
|
lambda x: np.exp(-x**2),
|
|
lambda x: np.ones_like(x)
|
|
]
|
|
|
|
def random_color(min_saturation=0.5, lightness_range=(0.5, 1.0)):
|
|
''' Returns a random color where the maximum lightness is 0.5 and the saturation is always greater than 0.5. '''
|
|
h = np.random.uniform(0, 360) # Hue is a value between 0 and 360
|
|
s = np.random.uniform(min_saturation, 1) # Saturation is a value between min_saturation and 1
|
|
l = np.random.uniform(lightness_range[0], lightness_range[1]) # Lightness is a value between 0 and max_lightness
|
|
return [x for x in colorsys.hls_to_rgb(h/360, l, s)]# Convert HSL to RGB
|
|
|
|
def test_cubic_lagrange_polynomial():
|
|
x = np.linspace(0, 16, 100)
|
|
x1, x2, x3, x4 = 1, 3, 7, 15
|
|
y1, y2, y3, y4 = 1, -2, 3.5, 5
|
|
|
|
a0, a1, a2, a3 = polyfit4(x1, x2, x3, x4, y1, y2, y3, y4)
|
|
print(a0, a1, a2, a3)
|
|
y = cubic_poly(x, a0, a1, a2, a3)
|
|
|
|
# plt.plot(x, cubic_lagrange(x, x1, x2, x3, x4, y1, y2, y3, y4), 'b', label='cubic Lagrange ChatGPT one shot')
|
|
plt.plot(x, y, 'g', label='Cubic Lagrange polynomial through linear least-squares')
|
|
plt.plot([x1, x2, x3, x4], [y1, y2, y3, y4], 'ro')
|
|
plt.grid()
|
|
plt.show()
|
|
|
|
def test_cubic_lagrange_minimize():
|
|
for f in functions_list:
|
|
a, b = -1.2, 1.5
|
|
x = np.linspace(a, b, 100)
|
|
x_min = cubic_lagrange_minimize(f, a, b, tol=1e-6, callback=cubic_lagrange_minimize_callback_detailed)
|
|
print('Solution : ', x_min, f(x_min))
|
|
c = random_color(.5, (0.5, 0.7))
|
|
plt.plot(x, f(x), color=c)
|
|
plt.plot(x_min, f(x_min), 'o', color=c)
|
|
plt.grid()
|
|
plt.show()
|
|
|
|
# test_cubic_lagrange_polynomial()
|
|
test_cubic_lagrange_minimize() |