Source code for PySpice.Math.Calculus

####################################################################################################
#
# PySpice - A Spice Package for Python
# Copyright (C) 2014 Fabrice Salvaire
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
####################################################################################################

"""This module provides algorithms to compute the derivative of a function sampled on an uniform
grid.
"""

####################################################################################################

import fractions

import numpy as np

####################################################################################################

from PySpice.Math import odd

####################################################################################################

[docs]def compute_exact_finite_difference_coefficients(derivative_order, grid, x0=0): """This function compute the finite difference coefficients for the given derivative order and grid. The parameter *x* specifies where is computed the derivative on the grid. The grid is given as a list of integer offsets. The algorithm is derived from the article: Generation of Finite Difference Formulas on Arbitrary Space Grids, Bengt Fornberg, Mathematics of computation, volume 51, number 184, october 1988 """ N = len(grid) # d[m,n,v] d = [[[0 for v in range(N)] for n in range(N)] for m in range(derivative_order +1)] d[0][0][0] = fractions.Fraction(1,1) c1 = 1 for n in range(1, N): c2 = 1 for v in range(n): c3 = grid[n] - grid[v] c2 *= c3 if n <= derivative_order: d[n][n-1][v] = 0 for m in range(min(n, derivative_order) +1): d[m][n][v] = ( (grid[n] - x0)*d[m][n-1][v] - m*d[m-1][n-1][v] ) / c3 for m in range(min(n, derivative_order) +1): d[m][n][n] = fractions.Fraction(c1,c2)*( m*d[m-1][n-1][n-1] - (grid[n-1] - x0)*d[m][n-1][n-1] ) c1 = c2 return d[-1][-1]
####################################################################################################
[docs]def compute_finite_difference_coefficients(derivative_order, grid): return [float(x) for x in compute_exact_finite_difference_coefficients(derivative_order, grid)]
#################################################################################################### _coefficient_cache = dict(centred={}, forward={}, backward={})
[docs]def get_finite_difference_coefficients(derivative_order, accuracy_order, grid_type): if derivative_order < 1: raise ValueError("Wrong derivative order") if odd(accuracy_order) or accuracy_order < 2: raise ValueError("Wrong accuracy order") if grid_type == 'centred': window_size = accuracy_order // 2 grid = list(range(-window_size, window_size +1)) elif grid_type == 'forward': grid = list(range(derivative_order + accuracy_order)) elif grid_type == 'backward': grid = list(range(-(derivative_order + accuracy_order) +1, 1)) grid = list(reversed(grid)) # Fixme: why ? else: raise ValueError("Wrong grid type") key = '{}-{}'.format(derivative_order, accuracy_order) coefficients = _coefficient_cache[grid_type].get(key, None) if coefficients is None: coefficients = compute_finite_difference_coefficients(derivative_order, grid) _coefficient_cache[grid_type][key] = coefficients return grid, coefficients
####################################################################################################
[docs]def simple_derivative(x, values): """ Compute the derivative as a simple slope. """ return x[:-1], np.diff(values)/np.diff(x)
####################################################################################################
[docs]def derivative(x, values, derivative_order=1, accuracy_order=4): """Compute the derivative at the given derivative order and accuracy order. The precision of the Taylor expansion is :math:`\mathcal{O}(dx^{accuracy})`. """ dx = np.diff(x) # if not np.all(dx == dx[0]): # raise ValueError("Sampling is not uniform") dx = dx[0] values_size, = values.shape derivative = np.zeros(values_size, dtype=values.dtype) grid, coefficients = get_finite_difference_coefficients(derivative_order, accuracy_order, 'centred') window_size = grid[-1] # print grid, coefficients vector_size = values_size - 2*window_size if not vector_size: raise ValueError("The size of the value's array is not sufficient for the given accuracy order") lower_index = window_size upper_index = values_size - window_size derivative_view = derivative[window_size:-window_size] for offset, coefficient in zip(grid, coefficients): if coefficient: # print offset, lower_index + offset, upper_index + offset derivative_view += values[lower_index + offset:upper_index + offset] * coefficient grid, coefficients = get_finite_difference_coefficients(derivative_order, accuracy_order, 'forward') # print grid, coefficients grid_size = len(grid) upper_index = window_size derivative_view = derivative[:window_size] for offset, coefficient in zip(grid, coefficients): # print offset, offset, window_size+offset derivative_view += values[offset:upper_index + offset] * coefficient grid, coefficients = get_finite_difference_coefficients(derivative_order, accuracy_order, 'backward') # print grid, coefficients grid_size = len(grid) lower_index = values_size - window_size upper_index = values_size derivative_view = derivative[-window_size:] for offset, coefficient in zip(grid, coefficients): # print offset, lower_index + offset, upper_index + offset derivative_view += values[lower_index + offset:upper_index + offset] * coefficient return derivative / dx**derivative_order