2
Finite difference weights
3
=========================
5
This module implements an algorithm for efficient generation of finite
6
difference weights for ordinary differentials of functions for
7
derivatives from 0 (interpolation) up to arbitrary order.
9
The core algorithm is provided in the finite difference weight generating
10
function (finite_diff_weights), and two convenience functions are provided
13
- estimating a derivative (or interpolate) directly from a series of points
14
is also provided (``apply_finite_diff``).
15
- making a finite difference approximation of a Derivative instance
21
from sympy.core.compatibility import iterable
24
def finite_diff_weights(order, x_list, x0):
26
Calculates the finite difference weights for an arbitrarily
27
spaced one-dimensional grid (x_list) for derivatives at 'x0'
28
of order 0, 1, ..., up to 'order' using a recursive formula.
34
Up to what derivative order weights should be calculated.
35
0 corresponds to interpolation.
37
Strictly monotonically increasing sequence of values for
38
the independent variable.
40
At what value of the independent variable the finite difference
41
weights should be generated.
47
A list of sublists, each corresponding to coefficients for
48
increasing derivative order, and each containing lists of
49
coefficients for increasing accuracy.
54
>>> from sympy import S
55
>>> from sympy.calculus import finite_diff_weights
56
>>> finite_diff_weights(1, [-S(1)/2, S(1)/2, S(3)/2, S(5)/2], 0)
60
[5/16, 15/16, -5/16, 1/16]],
61
[[0, 0, 0, 0], [-1, 1, 0, 0], [-1, 1, 0, 0], [-23/24, 7/8, 1/8, -1/24]]]
63
the result is two subslists, the first is for the 0:th derivative
64
(interpolation) and the second for the first derivative (we gave
65
1 as the parameter of order so this is why we get no list for
66
a higher order derivative). Each sublist contains the most accurate
67
formula in the end (all points used).
69
Beware of the offset in the lower accuracy formulae when looking at a
72
>>> from sympy import S
73
>>> from sympy.calculus import finite_diff_weights
74
>>> finite_diff_weights(1, [-S(5)/2, -S(3)/2, -S(1)/2, S(1)/2,
75
... S(3)/2, S(5)/2], 0) #doctest: +NORMALIZE_WHITESPACE
77
[-3/2, 5/2, 0, 0, 0, 0],
78
[3/8, -5/4, 15/8, 0, 0, 0],
79
[1/16, -5/16, 15/16, 5/16, 0, 0],
80
[3/128, -5/32, 45/64, 15/32, -5/128, 0],
81
[3/256, -25/256, 75/128, 75/128, -25/256, 3/256]],
85
[1/24, -1/8, -7/8, 23/24, 0, 0],
86
[0, 1/24, -9/8, 9/8, -1/24, 0],
87
[-3/640, 25/384, -75/64, 75/64, -25/384, 3/640]]]
90
The capability to generate weights at arbitrary points can be
91
used e.g. to minimize Runge's phenomenon by using Chebyshev nodes:
93
>>> from sympy import cos, symbols, pi, simplify
94
>>> from sympy.calculus import finite_diff_weights
95
>>> N, (h, x) = 4, symbols('h x')
96
>>> x_list = [x+h*cos(i*pi/(N)) for i in range(N,-1,-1)] # chebyshev nodes
98
[-h + x, -sqrt(2)*h/2 + x, x, sqrt(2)*h/2 + x, h + x]
99
>>> mycoeffs = finite_diff_weights(1, x_list, 0)[1][4]
100
>>> [simplify(c) for c in mycoeffs] #doctest: +NORMALIZE_WHITESPACE
101
[(h**3/2 + h**2*x - 3*h*x**2 - 4*x**3)/h**4,
102
(-sqrt(2)*h**3 - 4*h**2*x + 3*sqrt(2)*h*x**2 + 8*x**3)/h**4,
103
6*x/h**2 - 8*x**3/h**4,
104
(sqrt(2)*h**3 - 4*h**2*x - 3*sqrt(2)*h*x**2 + 8*x**3)/h**4,
105
(-h**3/2 + h**2*x + 3*h*x**2 - 4*x**3)/h**4]
110
If weights for a finite difference approximation
111
of the 3rd order derivative is wanted, weights for 0th, 1st
112
and 2nd order are calculated "for free", so are formulae using
113
fewer and fewer of the parameters. This is something one can
114
take advantage of to save computational cost.
119
sympy.calculus.finite_diff.apply_finite_diff
125
.. [1] Generation of Finite Difference Formulas on Arbitrarily Spaced
126
Grids, Bengt Fornberg; Mathematics of computation; 51; 184;
127
(1988); 699-706; doi:10.1090/S0025-5718-1988-0935077-0
130
# The notation below closely corresponds to the one used in the paper.
132
raise ValueError("Negative derivative order illegal.")
133
if int(order) != order:
134
raise ValueError("Non-integer order illegal")
137
delta = [[[0 for nu in range(N+1)] for n in range(N+1)] for
139
delta[0][0][0] = S(1)
141
for n in range(1, N+1):
143
for nu in range(0, n):
144
c3 = x_list[n]-x_list[nu]
147
delta[n][n-1][nu] = 0
148
for m in range(0, min(n, M)+1):
149
delta[m][n][nu] = (x_list[n]-x0)*delta[m][n-1][nu] -\
150
m*delta[m-1][n-1][nu]
151
delta[m][n][nu] /= c3
152
for m in range(0, min(n, M)+1):
153
delta[m][n][n] = c1/c2*(m*delta[m-1][n-1][n-1] -
154
(x_list[n-1]-x0)*delta[m][n-1][n-1])
159
def apply_finite_diff(order, x_list, y_list, x0):
161
Calculates the finite difference approximation of
162
the derivative of requested order at x0 from points
163
provided in x_list and y_list.
169
order of derivative to approximate. 0 corresponds to interpolation.
171
Strictly monotonically increasing sequence of values for
172
the independent variable.
174
The function value at corresponding values for the independent
177
At what value of the independent variable the derivative should be
183
sympy.core.add.Add or sympy.core.numbers.Number
184
The finite difference expression approximating the requested
185
derivative order at x0.
190
>>> from sympy.calculus import apply_finite_diff
191
>>> cube = lambda arg: (1.0*arg)**3
192
>>> xlist = range(-3,3+1)
193
>>> apply_finite_diff(2, xlist, map(cube, xlist), 2) - 12 # doctest: +SKIP
194
-3.55271367880050e-15
196
we see that the example above only contain rounding errors.
197
apply_finite_diff can also be used on more abstract objects:
199
>>> from sympy import IndexedBase, Idx
200
>>> from sympy.calculus import apply_finite_diff
201
>>> x, y = map(IndexedBase, 'xy')
203
>>> x_list, y_list = zip(*[(x[i+j], y[i+j]) for j in range(-1,2)])
204
>>> apply_finite_diff(1, x_list, y_list, x[i])
205
(-1 + (x[i + 1] - x[i])/(-x[i - 1] + x[i]))*y[i]/(x[i + 1] - x[i]) + \
206
(-x[i - 1] + x[i])*y[i + 1]/((-x[i - 1] + x[i + 1])*(x[i + 1] - x[i])) - \
207
(x[i + 1] - x[i])*y[i - 1]/((-x[i - 1] + x[i + 1])*(-x[i - 1] + x[i]))
213
Order = 0 corresponds to interpolation.
214
Only supply so many points you think makes sense
215
to around x0 when extracting the derivative (the function
216
need to be well behaved within that region). Also beware
217
of Runge's phenomenon.
222
sympy.calculus.finite_diff.finite_diff_weights
227
Fortran 90 implementation with Python interface for numerics: finitediff_
229
.. _finitediff: https://github.com/bjodah/finitediff
233
# In the original paper the following holds for the notation:
235
# N = len(x_list) - 1
238
if len(x_list) != len(y_list):
239
raise ValueError("x_list and y_list not equal in length.")
241
delta = finite_diff_weights(order, x_list, x0)
244
for nu in range(0, len(x_list)):
245
derivative += delta[order][N][nu]*y_list[nu]
249
def as_finite_diff(derivative, points=1, x0=None, wrt=None):
251
Returns an approximation of a derivative of a function in
252
the form of a finite difference formula. The expression is a
253
weighted sum of the function at a number of discrete values of
254
(one of) the independent variable(s).
259
derivative: a Derivative instance (needs to have an variables
262
points: sequence or coefficient, optional
263
If sequence: discrete values (length >= order+1) of the
264
independent variable used for generating the finite
266
If it is a coefficient, it will be used as the step-size
267
for generating an equidistant sequence of length order+1
268
centered around x0. default: 1 (step-size 1)
270
x0: number or Symbol, optional
271
the value of the independent variable (wrt) at which the
272
derivative is to be approximated. default: same as wrt
274
wrt: Symbol, optional
275
"with respect to" the variable for which the (partial)
276
derivative is to be approximated for. If not provided it
277
is required that the Derivative is ordinary. default: None
283
>>> from sympy import symbols, Function, exp, sqrt, Symbol, as_finite_diff
284
>>> x, h = symbols('x h')
285
>>> f = Function('f')
286
>>> as_finite_diff(f(x).diff(x))
287
-f(x - 1/2) + f(x + 1/2)
289
The default step size and number of points are 1 and ``order + 1``
290
respectively. We can change the step size by passing a symbol
293
>>> as_finite_diff(f(x).diff(x), h)
294
-f(-h/2 + x)/h + f(h/2 + x)/h
296
We can also specify the discretized values to be used in a sequence:
298
>>> as_finite_diff(f(x).diff(x), [x, x+h, x+2*h])
299
-3*f(x)/(2*h) + 2*f(h + x)/h - f(2*h + x)/(2*h)
301
The algorithm is not restricted to use equidistant spacing, nor
302
do we need to make the approximation around x0, but we can get
303
an expression estimating the derivative at an offset:
305
>>> e, sq2 = exp(1), sqrt(2)
306
>>> xl = [x-h, x+h, x+e*h]
307
>>> as_finite_diff(f(x).diff(x, 1), xl, x+h*sq2)
308
2*h*((h + sqrt(2)*h)/(2*h) - (-sqrt(2)*h + h)/(2*h))*f(E*h + x)/\
309
((-h + E*h)*(h + E*h)) + (-(-sqrt(2)*h + h)/(2*h) - \
310
(-sqrt(2)*h + E*h)/(2*h))*f(-h + x)/(h + E*h) + \
311
(-(h + sqrt(2)*h)/(2*h) + (-sqrt(2)*h + E*h)/(2*h))*f(h + x)/(-h + E*h)
313
Partial derivatives are also supported:
316
>>> d2fdxdy=f(x,y).diff(x,y)
317
>>> as_finite_diff(d2fdxdy, wrt=x)
318
-f(x - 1/2, y) + f(x + 1/2, y)
323
sympy.calculus.finite_diff.apply_finite_diff
324
sympy.calculus.finite_diff.finite_diff_weights
328
wrt = derivative.variables[0]
329
# we need Derivative to be univariate to guess wrt
330
if any(v != wrt for v in derivative.variables):
331
raise ValueError('if the function is not univariate' +
332
' then `wrt` must be given')
334
order = derivative.variables.count(wrt)
339
if not iterable(points):
340
# points is simply the step-size, let's make it a
341
# equidistant sequence centered around x0
343
# even order => odd number of points, grid point included
344
points = [x0 + points*i for i
345
in range(-order//2, order//2 + 1)]
347
# odd order => even number of points, half-way wrt grid point
348
points = [x0 + points*i/S(2) for i
349
in range(-order, order + 1, 2)]
351
if len(points) < order+1:
352
raise ValueError("Too few points for order %d" % order)
353
return apply_finite_diff(order, points, [
354
derivative.expr.subs({wrt: x}) for x in points], x0)
added
removed
lib/py/sympy/calculus/finite_diff.py
