Homec4science
Diffusion Public Numerical Analysis Jupyter Notebook (Modifications-for-2025) Browse

PubNumAnalysisIpynb/InterpolationLib.py
58eed21d4fd6Modifications-for-2025

Last Change
Raw File

InterpolationLib.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 30
@author: Simone Deparis
"""
import numpy as np
from types import FunctionType
# Defining the Lagrange basis functions
def LagrangeBasis(t,k,z):
# the input variables are:
# t : the interpolatory points
# k : which basis function
# z : where to evaluate the function
# careful, there are n+1 interpolation points!
n = t.size - 1
# init result to one, of same type and size as z
result = np.ones_like(z)
# first few checks on k:
if (type(k) != int) or (t.size < 1) or (k > n) or (k < 0):
raise ValueError('Lagrange basis needs a positive integer k, smaller than the size of t')
# loop on n to compute the product
for j in range(n+1) :
if (j == k) :
continue
if (t[k] == t[j]) :
raise ValueError('Lagrange basis: all the interpolation points need to be distinct')
result = result * (z - t[j]) / (t[k] - t[j])
return result
# Lagrange Interpolation of data (t,y), evaluated at ordinate(s) z
def LagrangeInterpolation(t,y,z):
# the input variables are:
# t : the interpolatory points
# y : the corresponding data at the points t
# z : where to evaluate the function
# {phi(t,k,.), k=0,...,n} is a basis of the polynomials of degree n
# y represents the coordinates of the interpolating polynomial with respect to this basis.
# Therefore LagrangePi(t,y,.) = y[0] phi(t,0,.) + ... + y[n] phi(t,n,.)
# careful, there are n+1 basis functions!
n = t.size - 1
# init result to zero, of same type and size as z
result = np.zeros_like(z)
# loop on n to compute the product
for k in range(0,n+1) :
result = result + y[k] * LagrangeBasis(t,k,z)
return result
# Piece-wise linear interpolation, equidistributed points
def PiecewiseLinearInterpolation(a,b,N,y,z):
# the input variables are:
# a,b : x[0] = a, x[n] = b
# N: the number of sub-intervals
# y : the corresponding data at the points x or the function f
# z : where to evaluate the function
def localLinearInterpolation (a, H, x, y, zk):
# if we are outside the interpolation interval, return the value at the closest node
if zk <= x[0]:
return x[0]
if zk >= x[-1]:
return x[-1]
# first find out in which interval we are. Easy here since equidistributed point
i = int( (zk-a)//H )
# if we are exactly at the node, return nodal value
if x[i] == zk :
return y[i]
# otherwise use formula for local linear interpolation
return y[i] + (y[i+1] - y[i])/(x[i+1] - x[i])* (zk - x[i])
# careful, there are N+1 points
H = (b-a)/N
if np.isclose(H,0) :
raise ValueError('PiecewiseLinearInterpolation : a and b are too close')
x = np.linspace(a,b,N+1)
# if y is a function, we need to evaluate it at all the locations in x
if isinstance(y, FunctionType) :
y = y(x)
# if we want to evaluate the interpolating function in a single point
if np.ndim(z) == 0 :
return localLinearInterpolation (b, H, x, y, z)
# if we are here, it means that we need to evaluate the interpolation on many points.
# init result to zero, of same type and size as z
result = np.zeros_like(z)
# The following part is not optimized
for k in range(z.size) :
result[k] = localLinearInterpolation (a, H, x, y, z[k])
return result
# Piece-wise quadratic interpolation, equidistributed points
def PiecewiseQuadraticInterpolation(a,b,N,f,z):
# the input variables are:
# a,b : x[0] = a, x[n] = b
# N: the number of sub-intervals
# f : the function to interpolate
# z : where to evaluate the function
assert isinstance(y, FunctionType), "f must be a function, not a list of values!"
def localQuadraticInterpolation (a, H, x, f, zk):
# if we are outside the interpolation interval, return the value at the closest node
if zk <= x[0]:
return x[0]
if zk >= x[-1]:
return x[-1]
# first find out in which interval we are. Easy here since equidistributed point
i = int( (zk-a)//H )
# if we are exactly at the node, return nodal value
if x[i] == zk :
return f(x[i])
t = np.array([ x[i], (x[i]+x[i+1])/2, x[i+1]] )
# otherwise use formula for local quadratic interpolation
return LagrangeInterpolation(t,f(t),zk)
# careful, there are N+1 points
H = (b-a)/N
if np.isclose(H,0) :
raise ValueError('PiecewiseQuadraticInterpolation : a and b are too close')
x = np.linspace(a,b,N+1)
# if we want to evaluate the interpolating function in a single point
if np.ndim(z) == 0 :
return localQuadraticInterpolation (b, H, x, f, z)
# if we are here, it means that we need to evaluate the interpolation on many points.
# init result to zero, of same type and size as z
result = np.zeros_like(z)
# The following part is not optimized
for k in range(z.size) :
result[k] = localQuadraticInterpolation (a, H, x, f, z[k])
return result
c4science · Help