图片解析应用
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

168 lines
7.6 KiB

  1. from itertools import product
  2. from sympy.core.function import (Function, diff)
  3. from sympy.core.numbers import Rational
  4. from sympy.core.singleton import S
  5. from sympy.core.symbol import symbols
  6. from sympy.functions.elementary.exponential import exp
  7. from sympy.calculus.finite_diff import (
  8. apply_finite_diff, differentiate_finite, finite_diff_weights,
  9. _as_finite_diff
  10. )
  11. from sympy.testing.pytest import raises, warns_deprecated_sympy
  12. def test_apply_finite_diff():
  13. x, h = symbols('x h')
  14. f = Function('f')
  15. assert (apply_finite_diff(1, [x-h, x+h], [f(x-h), f(x+h)], x) -
  16. (f(x+h)-f(x-h))/(2*h)).simplify() == 0
  17. assert (apply_finite_diff(1, [5, 6, 7], [f(5), f(6), f(7)], 5) -
  18. (Rational(-3, 2)*f(5) + 2*f(6) - S.Half*f(7))).simplify() == 0
  19. raises(ValueError, lambda: apply_finite_diff(1, [x, h], [f(x)]))
  20. def test_finite_diff_weights():
  21. d = finite_diff_weights(1, [5, 6, 7], 5)
  22. assert d[1][2] == [Rational(-3, 2), 2, Rational(-1, 2)]
  23. # Table 1, p. 702 in doi:10.1090/S0025-5718-1988-0935077-0
  24. # --------------------------------------------------------
  25. xl = [0, 1, -1, 2, -2, 3, -3, 4, -4]
  26. # d holds all coefficients
  27. d = finite_diff_weights(4, xl, S.Zero)
  28. # Zeroeth derivative
  29. for i in range(5):
  30. assert d[0][i] == [S.One] + [S.Zero]*8
  31. # First derivative
  32. assert d[1][0] == [S.Zero]*9
  33. assert d[1][2] == [S.Zero, S.Half, Rational(-1, 2)] + [S.Zero]*6
  34. assert d[1][4] == [S.Zero, Rational(2, 3), Rational(-2, 3), Rational(-1, 12), Rational(1, 12)] + [S.Zero]*4
  35. assert d[1][6] == [S.Zero, Rational(3, 4), Rational(-3, 4), Rational(-3, 20), Rational(3, 20),
  36. Rational(1, 60), Rational(-1, 60)] + [S.Zero]*2
  37. assert d[1][8] == [S.Zero, Rational(4, 5), Rational(-4, 5), Rational(-1, 5), Rational(1, 5),
  38. Rational(4, 105), Rational(-4, 105), Rational(-1, 280), Rational(1, 280)]
  39. # Second derivative
  40. for i in range(2):
  41. assert d[2][i] == [S.Zero]*9
  42. assert d[2][2] == [-S(2), S.One, S.One] + [S.Zero]*6
  43. assert d[2][4] == [Rational(-5, 2), Rational(4, 3), Rational(4, 3), Rational(-1, 12), Rational(-1, 12)] + [S.Zero]*4
  44. assert d[2][6] == [Rational(-49, 18), Rational(3, 2), Rational(3, 2), Rational(-3, 20), Rational(-3, 20),
  45. Rational(1, 90), Rational(1, 90)] + [S.Zero]*2
  46. assert d[2][8] == [Rational(-205, 72), Rational(8, 5), Rational(8, 5), Rational(-1, 5), Rational(-1, 5),
  47. Rational(8, 315), Rational(8, 315), Rational(-1, 560), Rational(-1, 560)]
  48. # Third derivative
  49. for i in range(3):
  50. assert d[3][i] == [S.Zero]*9
  51. assert d[3][4] == [S.Zero, -S.One, S.One, S.Half, Rational(-1, 2)] + [S.Zero]*4
  52. assert d[3][6] == [S.Zero, Rational(-13, 8), Rational(13, 8), S.One, -S.One,
  53. Rational(-1, 8), Rational(1, 8)] + [S.Zero]*2
  54. assert d[3][8] == [S.Zero, Rational(-61, 30), Rational(61, 30), Rational(169, 120), Rational(-169, 120),
  55. Rational(-3, 10), Rational(3, 10), Rational(7, 240), Rational(-7, 240)]
  56. # Fourth derivative
  57. for i in range(4):
  58. assert d[4][i] == [S.Zero]*9
  59. assert d[4][4] == [S(6), -S(4), -S(4), S.One, S.One] + [S.Zero]*4
  60. assert d[4][6] == [Rational(28, 3), Rational(-13, 2), Rational(-13, 2), S(2), S(2),
  61. Rational(-1, 6), Rational(-1, 6)] + [S.Zero]*2
  62. assert d[4][8] == [Rational(91, 8), Rational(-122, 15), Rational(-122, 15), Rational(169, 60), Rational(169, 60),
  63. Rational(-2, 5), Rational(-2, 5), Rational(7, 240), Rational(7, 240)]
  64. # Table 2, p. 703 in doi:10.1090/S0025-5718-1988-0935077-0
  65. # --------------------------------------------------------
  66. xl = [[j/S(2) for j in list(range(-i*2+1, 0, 2))+list(range(1, i*2+1, 2))]
  67. for i in range(1, 5)]
  68. # d holds all coefficients
  69. d = [finite_diff_weights({0: 1, 1: 2, 2: 4, 3: 4}[i], xl[i], 0) for
  70. i in range(4)]
  71. # Zeroth derivative
  72. assert d[0][0][1] == [S.Half, S.Half]
  73. assert d[1][0][3] == [Rational(-1, 16), Rational(9, 16), Rational(9, 16), Rational(-1, 16)]
  74. assert d[2][0][5] == [Rational(3, 256), Rational(-25, 256), Rational(75, 128), Rational(75, 128),
  75. Rational(-25, 256), Rational(3, 256)]
  76. assert d[3][0][7] == [Rational(-5, 2048), Rational(49, 2048), Rational(-245, 2048), Rational(1225, 2048),
  77. Rational(1225, 2048), Rational(-245, 2048), Rational(49, 2048), Rational(-5, 2048)]
  78. # First derivative
  79. assert d[0][1][1] == [-S.One, S.One]
  80. assert d[1][1][3] == [Rational(1, 24), Rational(-9, 8), Rational(9, 8), Rational(-1, 24)]
  81. assert d[2][1][5] == [Rational(-3, 640), Rational(25, 384), Rational(-75, 64),
  82. Rational(75, 64), Rational(-25, 384), Rational(3, 640)]
  83. assert d[3][1][7] == [Rational(5, 7168), Rational(-49, 5120),
  84. Rational(245, 3072), Rational(-1225, 1024),
  85. Rational(1225, 1024), Rational(-245, 3072),
  86. Rational(49, 5120), Rational(-5, 7168)]
  87. # Reasonably the rest of the table is also correct... (testing of that
  88. # deemed excessive at the moment)
  89. raises(ValueError, lambda: finite_diff_weights(-1, [1, 2]))
  90. raises(ValueError, lambda: finite_diff_weights(1.2, [1, 2]))
  91. x = symbols('x')
  92. raises(ValueError, lambda: finite_diff_weights(x, [1, 2]))
  93. def test_as_finite_diff():
  94. x = symbols('x')
  95. f = Function('f')
  96. dx = Function('dx')
  97. _as_finite_diff(f(x).diff(x), [x-2, x-1, x, x+1, x+2])
  98. # Use of undefined functions in ``points``
  99. df_true = -f(x+dx(x)/2-dx(x+dx(x)/2)/2) / dx(x+dx(x)/2) \
  100. + f(x+dx(x)/2+dx(x+dx(x)/2)/2) / dx(x+dx(x)/2)
  101. df_test = diff(f(x), x).as_finite_difference(points=dx(x), x0=x+dx(x)/2)
  102. assert (df_test - df_true).simplify() == 0
  103. def test_differentiate_finite():
  104. x, y, h = symbols('x y h')
  105. f = Function('f')
  106. with warns_deprecated_sympy():
  107. res0 = differentiate_finite(f(x, y) + exp(42), x, y, evaluate=True)
  108. xm, xp, ym, yp = [v + sign*S.Half for v, sign in product([x, y], [-1, 1])]
  109. ref0 = f(xm, ym) + f(xp, yp) - f(xm, yp) - f(xp, ym)
  110. assert (res0 - ref0).simplify() == 0
  111. g = Function('g')
  112. with warns_deprecated_sympy():
  113. res1 = differentiate_finite(f(x)*g(x) + 42, x, evaluate=True)
  114. ref1 = (-f(x - S.Half) + f(x + S.Half))*g(x) + \
  115. (-g(x - S.Half) + g(x + S.Half))*f(x)
  116. assert (res1 - ref1).simplify() == 0
  117. res2 = differentiate_finite(f(x) + x**3 + 42, x, points=[x-1, x+1])
  118. ref2 = (f(x + 1) + (x + 1)**3 - f(x - 1) - (x - 1)**3)/2
  119. assert (res2 - ref2).simplify() == 0
  120. raises(TypeError, lambda: differentiate_finite(f(x)*g(x), x,
  121. pints=[x-1, x+1]))
  122. res3 = differentiate_finite(f(x)*g(x).diff(x), x)
  123. ref3 = (-g(x) + g(x + 1))*f(x + S.Half) - (g(x) - g(x - 1))*f(x - S.Half)
  124. assert res3 == ref3
  125. res4 = differentiate_finite(f(x)*g(x).diff(x).diff(x), x)
  126. ref4 = -((g(x - Rational(3, 2)) - 2*g(x - S.Half) + g(x + S.Half))*f(x - S.Half)) \
  127. + (g(x - S.Half) - 2*g(x + S.Half) + g(x + Rational(3, 2)))*f(x + S.Half)
  128. assert res4 == ref4
  129. res5_expr = f(x).diff(x)*g(x).diff(x)
  130. res5 = differentiate_finite(res5_expr, points=[x-h, x, x+h])
  131. ref5 = (-2*f(x)/h + f(-h + x)/(2*h) + 3*f(h + x)/(2*h))*(-2*g(x)/h + g(-h + x)/(2*h) \
  132. + 3*g(h + x)/(2*h))/(2*h) - (2*f(x)/h - 3*f(-h + x)/(2*h) - \
  133. f(h + x)/(2*h))*(2*g(x)/h - 3*g(-h + x)/(2*h) - g(h + x)/(2*h))/(2*h)
  134. assert res5 == ref5
  135. res6 = res5.limit(h, 0).doit()
  136. ref6 = diff(res5_expr, x)
  137. assert res6 == ref6