Coverage for /usr/lib/python3/dist-packages/sympy/calculus/finite_diff.py: 12%
82 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Finite difference weights
3=========================
5This module implements an algorithm for efficient generation of finite
6difference weights for ordinary differentials of functions for
7derivatives from 0 (interpolation) up to arbitrary order.
9The core algorithm is provided in the finite difference weight generating
10function (``finite_diff_weights``), and two convenience functions are provided
11for:
13- estimating a derivative (or interpolate) directly from a series of points
14 is also provided (``apply_finite_diff``).
15- differentiating by using finite difference approximations
16 (``differentiate_finite``).
18"""
20from sympy.core.function import Derivative
21from sympy.core.singleton import S
22from sympy.core.function import Subs
23from sympy.core.traversal import preorder_traversal
24from sympy.utilities.exceptions import sympy_deprecation_warning
25from sympy.utilities.iterables import iterable
29def finite_diff_weights(order, x_list, x0=S.One):
30 """
31 Calculates the finite difference weights for an arbitrarily spaced
32 one-dimensional grid (``x_list``) for derivatives at ``x0`` of order
33 0, 1, ..., up to ``order`` using a recursive formula. Order of accuracy
34 is at least ``len(x_list) - order``, if ``x_list`` is defined correctly.
36 Parameters
37 ==========
39 order: int
40 Up to what derivative order weights should be calculated.
41 0 corresponds to interpolation.
42 x_list: sequence
43 Sequence of (unique) values for the independent variable.
44 It is useful (but not necessary) to order ``x_list`` from
45 nearest to furthest from ``x0``; see examples below.
46 x0: Number or Symbol
47 Root or value of the independent variable for which the finite
48 difference weights should be generated. Default is ``S.One``.
50 Returns
51 =======
53 list
54 A list of sublists, each corresponding to coefficients for
55 increasing derivative order, and each containing lists of
56 coefficients for increasing subsets of x_list.
58 Examples
59 ========
61 >>> from sympy import finite_diff_weights, S
62 >>> res = finite_diff_weights(1, [-S(1)/2, S(1)/2, S(3)/2, S(5)/2], 0)
63 >>> res
64 [[[1, 0, 0, 0],
65 [1/2, 1/2, 0, 0],
66 [3/8, 3/4, -1/8, 0],
67 [5/16, 15/16, -5/16, 1/16]],
68 [[0, 0, 0, 0],
69 [-1, 1, 0, 0],
70 [-1, 1, 0, 0],
71 [-23/24, 7/8, 1/8, -1/24]]]
72 >>> res[0][-1] # FD weights for 0th derivative, using full x_list
73 [5/16, 15/16, -5/16, 1/16]
74 >>> res[1][-1] # FD weights for 1st derivative
75 [-23/24, 7/8, 1/8, -1/24]
76 >>> res[1][-2] # FD weights for 1st derivative, using x_list[:-1]
77 [-1, 1, 0, 0]
78 >>> res[1][-1][0] # FD weight for 1st deriv. for x_list[0]
79 -23/24
80 >>> res[1][-1][1] # FD weight for 1st deriv. for x_list[1], etc.
81 7/8
83 Each sublist contains the most accurate formula at the end.
84 Note, that in the above example ``res[1][1]`` is the same as ``res[1][2]``.
85 Since res[1][2] has an order of accuracy of
86 ``len(x_list[:3]) - order = 3 - 1 = 2``, the same is true for ``res[1][1]``!
88 >>> res = finite_diff_weights(1, [S(0), S(1), -S(1), S(2), -S(2)], 0)[1]
89 >>> res
90 [[0, 0, 0, 0, 0],
91 [-1, 1, 0, 0, 0],
92 [0, 1/2, -1/2, 0, 0],
93 [-1/2, 1, -1/3, -1/6, 0],
94 [0, 2/3, -2/3, -1/12, 1/12]]
95 >>> res[0] # no approximation possible, using x_list[0] only
96 [0, 0, 0, 0, 0]
97 >>> res[1] # classic forward step approximation
98 [-1, 1, 0, 0, 0]
99 >>> res[2] # classic centered approximation
100 [0, 1/2, -1/2, 0, 0]
101 >>> res[3:] # higher order approximations
102 [[-1/2, 1, -1/3, -1/6, 0], [0, 2/3, -2/3, -1/12, 1/12]]
104 Let us compare this to a differently defined ``x_list``. Pay attention to
105 ``foo[i][k]`` corresponding to the gridpoint defined by ``x_list[k]``.
107 >>> foo = finite_diff_weights(1, [-S(2), -S(1), S(0), S(1), S(2)], 0)[1]
108 >>> foo
109 [[0, 0, 0, 0, 0],
110 [-1, 1, 0, 0, 0],
111 [1/2, -2, 3/2, 0, 0],
112 [1/6, -1, 1/2, 1/3, 0],
113 [1/12, -2/3, 0, 2/3, -1/12]]
114 >>> foo[1] # not the same and of lower accuracy as res[1]!
115 [-1, 1, 0, 0, 0]
116 >>> foo[2] # classic double backward step approximation
117 [1/2, -2, 3/2, 0, 0]
118 >>> foo[4] # the same as res[4]
119 [1/12, -2/3, 0, 2/3, -1/12]
121 Note that, unless you plan on using approximations based on subsets of
122 ``x_list``, the order of gridpoints does not matter.
124 The capability to generate weights at arbitrary points can be
125 used e.g. to minimize Runge's phenomenon by using Chebyshev nodes:
127 >>> from sympy import cos, symbols, pi, simplify
128 >>> N, (h, x) = 4, symbols('h x')
129 >>> x_list = [x+h*cos(i*pi/(N)) for i in range(N,-1,-1)] # chebyshev nodes
130 >>> print(x_list)
131 [-h + x, -sqrt(2)*h/2 + x, x, sqrt(2)*h/2 + x, h + x]
132 >>> mycoeffs = finite_diff_weights(1, x_list, 0)[1][4]
133 >>> [simplify(c) for c in mycoeffs] #doctest: +NORMALIZE_WHITESPACE
134 [(h**3/2 + h**2*x - 3*h*x**2 - 4*x**3)/h**4,
135 (-sqrt(2)*h**3 - 4*h**2*x + 3*sqrt(2)*h*x**2 + 8*x**3)/h**4,
136 (6*h**2*x - 8*x**3)/h**4,
137 (sqrt(2)*h**3 - 4*h**2*x - 3*sqrt(2)*h*x**2 + 8*x**3)/h**4,
138 (-h**3/2 + h**2*x + 3*h*x**2 - 4*x**3)/h**4]
140 Notes
141 =====
143 If weights for a finite difference approximation of 3rd order
144 derivative is wanted, weights for 0th, 1st and 2nd order are
145 calculated "for free", so are formulae using subsets of ``x_list``.
146 This is something one can take advantage of to save computational cost.
147 Be aware that one should define ``x_list`` from nearest to furthest from
148 ``x0``. If not, subsets of ``x_list`` will yield poorer approximations,
149 which might not grand an order of accuracy of ``len(x_list) - order``.
151 See also
152 ========
154 sympy.calculus.finite_diff.apply_finite_diff
156 References
157 ==========
159 .. [1] Generation of Finite Difference Formulas on Arbitrarily Spaced
160 Grids, Bengt Fornberg; Mathematics of computation; 51; 184;
161 (1988); 699-706; doi:10.1090/S0025-5718-1988-0935077-0
163 """
164 # The notation below closely corresponds to the one used in the paper.
165 order = S(order)
166 if not order.is_number:
167 raise ValueError("Cannot handle symbolic order.")
168 if order < 0:
169 raise ValueError("Negative derivative order illegal.")
170 if int(order) != order:
171 raise ValueError("Non-integer order illegal")
172 M = order
173 N = len(x_list) - 1
174 delta = [[[0 for nu in range(N+1)] for n in range(N+1)] for
175 m in range(M+1)]
176 delta[0][0][0] = S.One
177 c1 = S.One
178 for n in range(1, N+1):
179 c2 = S.One
180 for nu in range(n):
181 c3 = x_list[n] - x_list[nu]
182 c2 = c2 * c3
183 if n <= M:
184 delta[n][n-1][nu] = 0
185 for m in range(min(n, M)+1):
186 delta[m][n][nu] = (x_list[n]-x0)*delta[m][n-1][nu] -\
187 m*delta[m-1][n-1][nu]
188 delta[m][n][nu] /= c3
189 for m in range(min(n, M)+1):
190 delta[m][n][n] = c1/c2*(m*delta[m-1][n-1][n-1] -
191 (x_list[n-1]-x0)*delta[m][n-1][n-1])
192 c1 = c2
193 return delta
196def apply_finite_diff(order, x_list, y_list, x0=S.Zero):
197 """
198 Calculates the finite difference approximation of
199 the derivative of requested order at ``x0`` from points
200 provided in ``x_list`` and ``y_list``.
202 Parameters
203 ==========
205 order: int
206 order of derivative to approximate. 0 corresponds to interpolation.
207 x_list: sequence
208 Sequence of (unique) values for the independent variable.
209 y_list: sequence
210 The function value at corresponding values for the independent
211 variable in x_list.
212 x0: Number or Symbol
213 At what value of the independent variable the derivative should be
214 evaluated. Defaults to 0.
216 Returns
217 =======
219 sympy.core.add.Add or sympy.core.numbers.Number
220 The finite difference expression approximating the requested
221 derivative order at ``x0``.
223 Examples
224 ========
226 >>> from sympy import apply_finite_diff
227 >>> cube = lambda arg: (1.0*arg)**3
228 >>> xlist = range(-3,3+1)
229 >>> apply_finite_diff(2, xlist, map(cube, xlist), 2) - 12 # doctest: +SKIP
230 -3.55271367880050e-15
232 we see that the example above only contain rounding errors.
233 apply_finite_diff can also be used on more abstract objects:
235 >>> from sympy import IndexedBase, Idx
236 >>> x, y = map(IndexedBase, 'xy')
237 >>> i = Idx('i')
238 >>> x_list, y_list = zip(*[(x[i+j], y[i+j]) for j in range(-1,2)])
239 >>> apply_finite_diff(1, x_list, y_list, x[i])
240 ((x[i + 1] - x[i])/(-x[i - 1] + x[i]) - 1)*y[i]/(x[i + 1] - x[i]) -
241 (x[i + 1] - x[i])*y[i - 1]/((x[i + 1] - x[i - 1])*(-x[i - 1] + x[i])) +
242 (-x[i - 1] + x[i])*y[i + 1]/((x[i + 1] - x[i - 1])*(x[i + 1] - x[i]))
244 Notes
245 =====
247 Order = 0 corresponds to interpolation.
248 Only supply so many points you think makes sense
249 to around x0 when extracting the derivative (the function
250 need to be well behaved within that region). Also beware
251 of Runge's phenomenon.
253 See also
254 ========
256 sympy.calculus.finite_diff.finite_diff_weights
258 References
259 ==========
261 Fortran 90 implementation with Python interface for numerics: finitediff_
263 .. _finitediff: https://github.com/bjodah/finitediff
265 """
267 # In the original paper the following holds for the notation:
268 # M = order
269 # N = len(x_list) - 1
271 N = len(x_list) - 1
272 if len(x_list) != len(y_list):
273 raise ValueError("x_list and y_list not equal in length.")
275 delta = finite_diff_weights(order, x_list, x0)
277 derivative = 0
278 for nu in range(len(x_list)):
279 derivative += delta[order][N][nu]*y_list[nu]
280 return derivative
283def _as_finite_diff(derivative, points=1, x0=None, wrt=None):
284 """
285 Returns an approximation of a derivative of a function in
286 the form of a finite difference formula. The expression is a
287 weighted sum of the function at a number of discrete values of
288 (one of) the independent variable(s).
290 Parameters
291 ==========
293 derivative: a Derivative instance
295 points: sequence or coefficient, optional
296 If sequence: discrete values (length >= order+1) of the
297 independent variable used for generating the finite
298 difference weights.
299 If it is a coefficient, it will be used as the step-size
300 for generating an equidistant sequence of length order+1
301 centered around ``x0``. default: 1 (step-size 1)
303 x0: number or Symbol, optional
304 the value of the independent variable (``wrt``) at which the
305 derivative is to be approximated. Default: same as ``wrt``.
307 wrt: Symbol, optional
308 "with respect to" the variable for which the (partial)
309 derivative is to be approximated for. If not provided it
310 is required that the Derivative is ordinary. Default: ``None``.
312 Examples
313 ========
315 >>> from sympy import symbols, Function, exp, sqrt, Symbol
316 >>> from sympy.calculus.finite_diff import _as_finite_diff
317 >>> x, h = symbols('x h')
318 >>> f = Function('f')
319 >>> _as_finite_diff(f(x).diff(x))
320 -f(x - 1/2) + f(x + 1/2)
322 The default step size and number of points are 1 and ``order + 1``
323 respectively. We can change the step size by passing a symbol
324 as a parameter:
326 >>> _as_finite_diff(f(x).diff(x), h)
327 -f(-h/2 + x)/h + f(h/2 + x)/h
329 We can also specify the discretized values to be used in a sequence:
331 >>> _as_finite_diff(f(x).diff(x), [x, x+h, x+2*h])
332 -3*f(x)/(2*h) + 2*f(h + x)/h - f(2*h + x)/(2*h)
334 The algorithm is not restricted to use equidistant spacing, nor
335 do we need to make the approximation around ``x0``, but we can get
336 an expression estimating the derivative at an offset:
338 >>> e, sq2 = exp(1), sqrt(2)
339 >>> xl = [x-h, x+h, x+e*h]
340 >>> _as_finite_diff(f(x).diff(x, 1), xl, x+h*sq2)
341 2*h*((h + sqrt(2)*h)/(2*h) - (-sqrt(2)*h + h)/(2*h))*f(E*h + x)/((-h + E*h)*(h + E*h)) +
342 (-(-sqrt(2)*h + h)/(2*h) - (-sqrt(2)*h + E*h)/(2*h))*f(-h + x)/(h + E*h) +
343 (-(h + sqrt(2)*h)/(2*h) + (-sqrt(2)*h + E*h)/(2*h))*f(h + x)/(-h + E*h)
345 Partial derivatives are also supported:
347 >>> y = Symbol('y')
348 >>> d2fdxdy=f(x,y).diff(x,y)
349 >>> _as_finite_diff(d2fdxdy, wrt=x)
350 -Derivative(f(x - 1/2, y), y) + Derivative(f(x + 1/2, y), y)
352 See also
353 ========
355 sympy.calculus.finite_diff.apply_finite_diff
356 sympy.calculus.finite_diff.finite_diff_weights
358 """
359 if derivative.is_Derivative:
360 pass
361 elif derivative.is_Atom:
362 return derivative
363 else:
364 return derivative.fromiter(
365 [_as_finite_diff(ar, points, x0, wrt) for ar
366 in derivative.args], **derivative.assumptions0)
368 if wrt is None:
369 old = None
370 for v in derivative.variables:
371 if old is v:
372 continue
373 derivative = _as_finite_diff(derivative, points, x0, v)
374 old = v
375 return derivative
377 order = derivative.variables.count(wrt)
379 if x0 is None:
380 x0 = wrt
382 if not iterable(points):
383 if getattr(points, 'is_Function', False) and wrt in points.args:
384 points = points.subs(wrt, x0)
385 # points is simply the step-size, let's make it a
386 # equidistant sequence centered around x0
387 if order % 2 == 0:
388 # even order => odd number of points, grid point included
389 points = [x0 + points*i for i
390 in range(-order//2, order//2 + 1)]
391 else:
392 # odd order => even number of points, half-way wrt grid point
393 points = [x0 + points*S(i)/2 for i
394 in range(-order, order + 1, 2)]
395 others = [wrt, 0]
396 for v in set(derivative.variables):
397 if v == wrt:
398 continue
399 others += [v, derivative.variables.count(v)]
400 if len(points) < order+1:
401 raise ValueError("Too few points for order %d" % order)
402 return apply_finite_diff(order, points, [
403 Derivative(derivative.expr.subs({wrt: x}), *others) for
404 x in points], x0)
407def differentiate_finite(expr, *symbols,
408 points=1, x0=None, wrt=None, evaluate=False):
409 r""" Differentiate expr and replace Derivatives with finite differences.
411 Parameters
412 ==========
414 expr : expression
415 \*symbols : differentiate with respect to symbols
416 points: sequence, coefficient or undefined function, optional
417 see ``Derivative.as_finite_difference``
418 x0: number or Symbol, optional
419 see ``Derivative.as_finite_difference``
420 wrt: Symbol, optional
421 see ``Derivative.as_finite_difference``
423 Examples
424 ========
426 >>> from sympy import sin, Function, differentiate_finite
427 >>> from sympy.abc import x, y, h
428 >>> f, g = Function('f'), Function('g')
429 >>> differentiate_finite(f(x)*g(x), x, points=[x-h, x+h])
430 -f(-h + x)*g(-h + x)/(2*h) + f(h + x)*g(h + x)/(2*h)
432 ``differentiate_finite`` works on any expression, including the expressions
433 with embedded derivatives:
435 >>> differentiate_finite(f(x) + sin(x), x, 2)
436 -2*f(x) + f(x - 1) + f(x + 1) - 2*sin(x) + sin(x - 1) + sin(x + 1)
437 >>> differentiate_finite(f(x, y), x, y)
438 f(x - 1/2, y - 1/2) - f(x - 1/2, y + 1/2) - f(x + 1/2, y - 1/2) + f(x + 1/2, y + 1/2)
439 >>> differentiate_finite(f(x)*g(x).diff(x), x)
440 (-g(x) + g(x + 1))*f(x + 1/2) - (g(x) - g(x - 1))*f(x - 1/2)
442 To make finite difference with non-constant discretization step use
443 undefined functions:
445 >>> dx = Function('dx')
446 >>> differentiate_finite(f(x)*g(x).diff(x), points=dx(x))
447 -(-g(x - dx(x)/2 - dx(x - dx(x)/2)/2)/dx(x - dx(x)/2) +
448 g(x - dx(x)/2 + dx(x - dx(x)/2)/2)/dx(x - dx(x)/2))*f(x - dx(x)/2)/dx(x) +
449 (-g(x + dx(x)/2 - dx(x + dx(x)/2)/2)/dx(x + dx(x)/2) +
450 g(x + dx(x)/2 + dx(x + dx(x)/2)/2)/dx(x + dx(x)/2))*f(x + dx(x)/2)/dx(x)
452 """
453 if any(term.is_Derivative for term in list(preorder_traversal(expr))):
454 evaluate = False
456 Dexpr = expr.diff(*symbols, evaluate=evaluate)
457 if evaluate:
458 sympy_deprecation_warning("""
459 The evaluate flag to differentiate_finite() is deprecated.
461 evaluate=True expands the intermediate derivatives before computing
462 differences, but this usually not what you want, as it does not
463 satisfy the product rule.
464 """,
465 deprecated_since_version="1.5",
466 active_deprecations_target="deprecated-differentiate_finite-evaluate",
467 )
468 return Dexpr.replace(
469 lambda arg: arg.is_Derivative,
470 lambda arg: arg.as_finite_difference(points=points, x0=x0, wrt=wrt))
471 else:
472 DFexpr = Dexpr.as_finite_difference(points=points, x0=x0, wrt=wrt)
473 return DFexpr.replace(
474 lambda arg: isinstance(arg, Subs),
475 lambda arg: arg.expr.as_finite_difference(
476 points=points, x0=arg.point[0], wrt=arg.variables[0]))