Calculus

Some calculus-related methods waiting to find a better place in the Diofant modules tree.

diofant.calculus.euler.euler_equations(L, funcs=(), vars=())[source]

Find the Euler-Lagrange equations for a given Lagrangian.

Parameters

  • L (Expr) – The Lagrangian that should be a function of the functions listed in the second argument and their derivatives.

    For example, in the case of two functions \(f(x,y)\), \(g(x,y)\) and two independent variables \(x\), \(y\) the Lagrangian would have the form:

    \[L\left(f(x,y),g(x,y),\frac{\partial f(x,y)}{\partial x}, \frac{\partial f(x,y)}{\partial y}, \frac{\partial g(x,y)}{\partial x}, \frac{\partial g(x,y)}{\partial y},x,y\right)\]

    In many cases it is not necessary to provide anything, except the Lagrangian, it will be auto-detected (and an error raised if this couldn’t be done).

  • funcs (Function or an iterable of Functions) – The functions that the Lagrangian depends on. The Euler equations are differential equations for each of these functions.

  • vars (Symbol or an iterable of Symbols) – The Symbols that are the independent variables of the functions.

Returns

eqns (list of Eq) – The list of differential equations, one for each function.

Examples

>>> L = (f(t).diff(t))**2/2 - f(t)**2/2
>>> euler_equations(L, f(t), t)
[Eq(-f(t) - Derivative(f(t), t, t), 0)]
>>> L = (f(t, x).diff(t))**2/2 - (f(t, x).diff(x))**2/2
>>> euler_equations(L, f(t, x), [t, x])
[Eq(-Derivative(f(t, x), t, t) + Derivative(f(t, x), x, x), 0)]

References

diofant.calculus.singularities.singularities(f, x)[source]

Find singularities of real-valued function \(f\) with respect to \(x\).

Examples

>>> singularities(1/(1 + x), x)
{-1}

>>> singularities(exp(1/x) + log(x + 1), x)
{-1, 0}

>>> singularities(exp(1/log(x + 1)), x)
{0}

Notes

Removable singularities are not supported now.

References

diofant.calculus.optimization.maximize(f, *v)[source]

Maximizes \(f\) with respect to given variables \(v\).

See also

minimize

diofant.calculus.optimization.minimize(f, *v)[source]

Minimizes \(f\) with respect to given variables \(v\).

Examples

>>> minimize(x**2, x)
(0, {x: 0})

>>> minimize([x**2, x >= 1], x)
(1, {x: 1})
>>> minimize([-x**2, x >= -2, x <= 1], x)
(-4, {x: -2})

See also

maximize

Finite difference weights

This module implements an algorithm for efficient generation of finite difference weights for ordinary differentials of functions for derivatives from 0 (interpolation) up to arbitrary order.

The core algorithm is provided in the finite difference weight generating function (finite_diff_weights), and two convenience functions are provided for:

  • estimating a derivative (or interpolate) directly from a series of points

    is also provided (apply_finite_diff).

  • making a finite difference approximation of a Derivative instance

    (as_finite_diff).

diofant.calculus.finite_diff.apply_finite_diff(order, x_list, y_list, x0=Integer(0))[source]

Calculates the finite difference approximation of the derivative of requested order at x0 from points provided in x_list and y_list.

Parameters

  • order (int) – order of derivative to approximate. 0 corresponds to interpolation.

  • x_list (sequence) – Sequence of (unique) values for the independent variable.

  • y_list (sequence) – The function value at corresponding values for the independent variable in x_list.

  • x0 (Number or Symbol) – At what value of the independent variable the derivative should be evaluated. Defaults to Integer(0).

Returns

diofant.core.add.Add or diofant.core.numbers.Number – The finite difference expression approximating the requested derivative order at x0.

Examples

>>> def cube(arg):
...     return (1.0*arg)**3
>>> xlist = range(-3, 4)
>>> apply_finite_diff(2, xlist, list(map(cube, xlist)), 2) - 12
-3.55271367880050e-15

we see that the example above only contain rounding errors. apply_finite_diff can also be used on more abstract objects:

>>> x, y = map(IndexedBase, 'xy')
>>> i = Idx('i')
>>> x_list, y_list = zip(*[(x[i + j], y[i + j]) for j in range(-1, 2)])
>>> apply_finite_diff(1, x_list, y_list, x[i])
(-1 + (x[i + 1] - x[i])/(-x[i - 1] + x[i]))*y[i]/(x[i + 1] - x[i]) +
(-x[i - 1] + x[i])*y[i + 1]/((-x[i - 1] + x[i + 1])*(x[i + 1] - x[i])) -
(x[i + 1] - x[i])*y[i - 1]/((-x[i - 1] + x[i + 1])*(-x[i - 1] + x[i]))

Notes

Order = 0 corresponds to interpolation. Only supply so many points you think makes sense to around x0 when extracting the derivative (the function need to be well behaved within that region). Also beware of Runge’s phenomenon.

See also

diofant.calculus.finite_diff.finite_diff_weights

References

Fortran 90 implementation with Python interface for numerics: finitediff

diofant.calculus.finite_diff.as_finite_diff(derivative, points=1, x0=None, wrt=None)[source]

Returns an approximation of a derivative of a function in the form of a finite difference formula. The expression is a weighted sum of the function at a number of discrete values of (one of) the independent variable(s).

Parameters

  • derivative (a Derivative instance (needs to have an variables) – and expr attribute).

  • points (sequence or coefficient, optional) – If sequence: discrete values (length >= order+1) of the independent variable used for generating the finite difference weights. If it is a coefficient, it will be used as the step-size for generating an equidistant sequence of length order+1 centered around x0. default: 1 (step-size 1)

  • x0 (number or Symbol, optional) – the value of the independent variable (wrt) at which the derivative is to be approximated. default: same as wrt

  • wrt (Symbol, optional) – “with respect to” the variable for which the (partial) derivative is to be approximated for. If not provided it is required that the Derivative is ordinary. default: None

Examples

>>> h = symbols('h')
>>> as_finite_diff(f(x).diff(x))
-f(x - 1/2) + f(x + 1/2)

The default step size and number of points are 1 and order + 1 respectively. We can change the step size by passing a symbol as a parameter:

>>> as_finite_diff(f(x).diff(x), h)
-f(-h/2 + x)/h + f(h/2 + x)/h

We can also specify the discretized values to be used in a sequence:

>>> as_finite_diff(f(x).diff(x), [x, x+h, x+2*h])
-3*f(x)/(2*h) + 2*f(h + x)/h - f(2*h + x)/(2*h)

The algorithm is not restricted to use equidistant spacing, nor do we need to make the approximation around x0, but we can get an expression estimating the derivative at an offset:

>>> e, sq2 = exp(1), sqrt(2)
>>> xl = [x-h, x+h, x+e*h]
>>> as_finite_diff(f(x).diff((x, 1)), xl, x+h*sq2)
2*h*f(E*h + x)*((h + sqrt(2)*h)/(2*h) -
(-sqrt(2)*h + h)/(2*h))/((-h + E*h)*(h + E*h)) +
f(-h + x)*(-(-sqrt(2)*h + h)/(2*h) - (-sqrt(2)*h + E*h)/(2*h))/(h +
E*h) + f(h + x)*(-(h + sqrt(2)*h)/(2*h) + (-sqrt(2)*h +
E*h)/(2*h))/(-h + E*h)

Partial derivatives are also supported:

>>> d2fdxdy = f(x, y).diff(x, y)
>>> as_finite_diff(d2fdxdy, wrt=x)
-f(x - 1/2, y) + f(x + 1/2, y)

See also

diofant.calculus.finite_diff.apply_finite_diff, diofant.calculus.finite_diff.finite_diff_weights

diofant.calculus.finite_diff.finite_diff_weights(order, x_list, x0=Integer(0))[source]

Calculates the finite difference weights for an arbitrarily spaced one-dimensional grid (x_list) for derivatives at ‘x0’ of order 0, 1, …, up to ‘order’ using a recursive formula. Order of accuracy is at least len(x_list) - order, if x_list is defined accurately.

Parameters

  • order (int) – Up to what derivative order weights should be calculated. 0 corresponds to interpolation.

  • x_list (sequence) – Sequence of (unique) values for the independent variable. It is useful (but not necessary) to order x_list from nearest to furthest from x0; see examples below.

  • x0 (Number or Symbol) – Root or value of the independent variable for which the finite difference weights should be generated. Defaults to Integer(0).

Returns

list – A list of sublists, each corresponding to coefficients for increasing derivative order, and each containing lists of coefficients for increasing subsets of x_list.

Examples

>>> res = finite_diff_weights(1, [-Rational(1, 2), Rational(1, 2),
...                               Rational(3, 2), Rational(5, 2)], 0)
>>> res
[[[1, 0, 0, 0],
  [1/2, 1/2, 0, 0],
  [3/8, 3/4, -1/8, 0],
  [5/16, 15/16, -5/16, 1/16]],
 [[0, 0, 0, 0],
  [-1, 1, 0, 0],
  [-1, 1, 0, 0],
  [-23/24, 7/8, 1/8, -1/24]]]
>>> res[0][-1]  # FD weights for 0th derivative, using full x_list
[5/16, 15/16, -5/16, 1/16]
>>> res[1][-1]  # FD weights for 1st derivative
[-23/24, 7/8, 1/8, -1/24]
>>> res[1][-2]  # FD weights for 1st derivative, using x_list[:-1]
[-1, 1, 0, 0]
>>> res[1][-1][0]  # FD weight for 1st deriv. for x_list[0]
-23/24
>>> res[1][-1][1]  # FD weight for 1st deriv. for x_list[1], etc.
7/8

Each sublist contains the most accurate formula at the end. Note, that in the above example res[1][1] is the same as res[1][2]. Since res[1][2] has an order of accuracy of len(x_list[:3]) - order = 3 - 1 = 2, the same is true for res[1][1]!

>>> res = finite_diff_weights(1, [Integer(0), Integer(1), -Integer(1),
...                               Integer(2), -Integer(2)], 0)[1]
>>> res
[[0, 0, 0, 0, 0],
 [-1, 1, 0, 0, 0],
 [0, 1/2, -1/2, 0, 0],
 [-1/2, 1, -1/3, -1/6, 0],
 [0, 2/3, -2/3, -1/12, 1/12]]
>>> res[0]  # no approximation possible, using x_list[0] only
[0, 0, 0, 0, 0]
>>> res[1]  # classic forward step approximation
[-1, 1, 0, 0, 0]
>>> res[2]  # classic centered approximation
[0, 1/2, -1/2, 0, 0]
>>> res[3:]  # higher order approximations
[[-1/2, 1, -1/3, -1/6, 0], [0, 2/3, -2/3, -1/12, 1/12]]

Let us compare this to a differently defined x_list. Pay attention to foo[i][k] corresponding to the gridpoint defined by x_list[k].

>>> foo = finite_diff_weights(1, [-Integer(2), -Integer(1), Integer(0),
...                               Integer(1), Integer(2)], 0)[1]
>>> foo
[[0, 0, 0, 0, 0],
 [-1, 1, 0, 0, 0],
 [1/2, -2, 3/2, 0, 0],
 [1/6, -1, 1/2, 1/3, 0],
 [1/12, -2/3, 0, 2/3, -1/12]]
>>> foo[1]  # not the same and of lower accuracy as res[1]!
[-1, 1, 0, 0, 0]
>>> foo[2]  # classic double backward step approximation
[1/2, -2, 3/2, 0, 0]
>>> foo[4]  # the same as res[4]
[1/12, -2/3, 0, 2/3, -1/12]

Note that, unless you plan on using approximations based on subsets of x_list, the order of gridpoints does not matter.

The capability to generate weights at arbitrary points can be used e.g. to minimize Runge’s phenomenon by using Chebyshev nodes:

>>> N, h = 4, symbols('h')
>>> x_list = [x + h*cos(i*pi/(N)) for i in range(N, -1, -1)]  # chebyshev nodes
>>> x_list
[-h + x, -sqrt(2)*h/2 + x, x, sqrt(2)*h/2 + x, h + x]
>>> mycoeffs = finite_diff_weights(1, x_list, 0)[1][4]
>>> [simplify(c) for c in mycoeffs]
[(h**3/2 + h**2*x - 3*h*x**2 - 4*x**3)/h**4,
(-sqrt(2)*h**3 - 4*h**2*x + 3*sqrt(2)*h*x**2 + 8*x**3)/h**4,
6*x/h**2 - 8*x**3/h**4,
(sqrt(2)*h**3 - 4*h**2*x - 3*sqrt(2)*h*x**2 + 8*x**3)/h**4,
(-h**3/2 + h**2*x + 3*h*x**2 - 4*x**3)/h**4]

Notes

If weights for a finite difference approximation of 3rd order derivative is wanted, weights for 0th, 1st and 2nd order are calculated “for free”, so are formulae using subsets of x_list. This is something one can take advantage of to save computational cost. Be aware that one should define x_list from nearest to farest from x_list. If not, subsets of x_list will yield poorer approximations, which might not grand an order of accuracy of len(x_list) - order.

See also

diofant.calculus.finite_diff.apply_finite_diff

References

  • Generation of Finite Difference Formulas on Arbitrarily Spaced Grids, Bengt Fornberg; Mathematics of computation; 51; 184; (1988); 699-706; doi:10.1090/S0025-5718-1988-0935077-0