Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85108220
InterpolationLib.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Sep 26, 20:18
Size
5 KB
Mime Type
text/x-python
Expires
Sat, Sep 28, 20:18 (2 d)
Engine
blob
Format
Raw Data
Handle
20896429
Attached To
rPUBNUMANALYSIS Public Numerical Analysis Jupyter Notebook
InterpolationLib.py
View Options
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 30
@author: Simone Deparis
"""
import
numpy
as
np
# 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
.
zeros_like
(
z
)
+
1
# 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
(
0
,
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
,
f
,
z
):
# the input variables are:
# a,b : x[0] = a, x[n] = b
# f : the corresponding data at the points x
# z : where to evaluate the function
def
localLinearInterpolation
(
a
,
H
,
x
,
y
,
zk
):
# first find out in which interval we are. Easy here since equidistributed point
i
=
int
(
(
zk
-
a
)
//
H
)
# if we are not in the interval, return constant function
if
x
[
i
]
==
zk
:
return
y
[
i
]
# 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 f is a function, we need to evaluate it at alle the locations
from
types
import
FunctionType
if
isinstance
(
f
,
FunctionType
)
:
y
=
f
(
x
)
else
:
y
=
f
# 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
# f : the corresponding data at the points x
# only works if fis a fonction
# z : where to evaluate the function
def
localQuadraticInterpolation
(
a
,
H
,
x
,
f
,
zk
):
# first find out in which interval we are. Easy here since equidistributed point
i
=
int
(
(
zk
-
a
)
//
H
)
# if we are not in the interval, return constant function
if
x
[
i
]
==
zk
:
return
f
(
x
[
i
])
t
=
np
.
array
([
x
[
i
],
(
x
[
i
]
+
x
[
i
+
1
])
/
2
,
x
[
i
+
1
]]
)
# 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 f is a function, we need to evaluate it at alle the locations
# from types import FunctionType
# if ~isinstance(f, FunctionType) :
# raise ValueError('PiecewiseQuadraticInterpolation : works only with a given function')
# 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
Event Timeline
Log In to Comment