test_pade.py 3.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107
  1. import numpy as np
  2. from scipy.interpolate import pade
  3. from scipy._lib._array_api import (
  4. xp_assert_equal, assert_array_almost_equal
  5. )
  6. def test_pade_trivial():
  7. nump, denomp = pade([1.0], 0)
  8. xp_assert_equal(nump.c, np.asarray([1.0]))
  9. xp_assert_equal(denomp.c, np.asarray([1.0]))
  10. nump, denomp = pade([1.0], 0, 0)
  11. xp_assert_equal(nump.c, np.asarray([1.0]))
  12. xp_assert_equal(denomp.c, np.asarray([1.0]))
  13. def test_pade_4term_exp():
  14. # First four Taylor coefficients of exp(x).
  15. # Unlike poly1d, the first array element is the zero-order term.
  16. an = [1.0, 1.0, 0.5, 1.0/6]
  17. nump, denomp = pade(an, 0)
  18. assert_array_almost_equal(nump.c, [1.0/6, 0.5, 1.0, 1.0])
  19. assert_array_almost_equal(denomp.c, [1.0])
  20. nump, denomp = pade(an, 1)
  21. assert_array_almost_equal(nump.c, [1.0/6, 2.0/3, 1.0])
  22. assert_array_almost_equal(denomp.c, [-1.0/3, 1.0])
  23. nump, denomp = pade(an, 2)
  24. assert_array_almost_equal(nump.c, [1.0/3, 1.0])
  25. assert_array_almost_equal(denomp.c, [1.0/6, -2.0/3, 1.0])
  26. nump, denomp = pade(an, 3)
  27. assert_array_almost_equal(nump.c, [1.0])
  28. assert_array_almost_equal(denomp.c, [-1.0/6, 0.5, -1.0, 1.0])
  29. # Testing inclusion of optional parameter
  30. nump, denomp = pade(an, 0, 3)
  31. assert_array_almost_equal(nump.c, [1.0/6, 0.5, 1.0, 1.0])
  32. assert_array_almost_equal(denomp.c, [1.0])
  33. nump, denomp = pade(an, 1, 2)
  34. assert_array_almost_equal(nump.c, [1.0/6, 2.0/3, 1.0])
  35. assert_array_almost_equal(denomp.c, [-1.0/3, 1.0])
  36. nump, denomp = pade(an, 2, 1)
  37. assert_array_almost_equal(nump.c, [1.0/3, 1.0])
  38. assert_array_almost_equal(denomp.c, [1.0/6, -2.0/3, 1.0])
  39. nump, denomp = pade(an, 3, 0)
  40. assert_array_almost_equal(nump.c, [1.0])
  41. assert_array_almost_equal(denomp.c, [-1.0/6, 0.5, -1.0, 1.0])
  42. # Testing reducing array.
  43. nump, denomp = pade(an, 0, 2)
  44. assert_array_almost_equal(nump.c, [0.5, 1.0, 1.0])
  45. assert_array_almost_equal(denomp.c, [1.0])
  46. nump, denomp = pade(an, 1, 1)
  47. assert_array_almost_equal(nump.c, [1.0/2, 1.0])
  48. assert_array_almost_equal(denomp.c, [-1.0/2, 1.0])
  49. nump, denomp = pade(an, 2, 0)
  50. assert_array_almost_equal(nump.c, [1.0])
  51. assert_array_almost_equal(denomp.c, [1.0/2, -1.0, 1.0])
  52. def test_pade_ints():
  53. # Simple test sequences (one of ints, one of floats).
  54. an_int = [1, 2, 3, 4]
  55. an_flt = [1.0, 2.0, 3.0, 4.0]
  56. # Make sure integer arrays give the same result as float arrays with same values.
  57. for i in range(0, len(an_int)):
  58. for j in range(0, len(an_int) - i):
  59. # Create float and int pade approximation for given order.
  60. nump_int, denomp_int = pade(an_int, i, j)
  61. nump_flt, denomp_flt = pade(an_flt, i, j)
  62. # Check that they are the same.
  63. xp_assert_equal(nump_int.c, nump_flt.c)
  64. xp_assert_equal(denomp_int.c, denomp_flt.c)
  65. def test_pade_complex():
  66. # Test sequence with known solutions - see page 6 of 10.1109/PESGM.2012.6344759.
  67. # Variable x is parameter - these tests will work with any complex number.
  68. x = 0.2 + 0.6j
  69. an = [1.0, x, -x*x.conjugate(), x.conjugate()*(x**2) + x*(x.conjugate()**2),
  70. -(x**3)*x.conjugate() - 3*(x*x.conjugate())**2 - x*(x.conjugate()**3)]
  71. nump, denomp = pade(an, 1, 1)
  72. assert_array_almost_equal(nump.c, [x + x.conjugate(), 1.0])
  73. assert_array_almost_equal(denomp.c, [x.conjugate(), 1.0])
  74. nump, denomp = pade(an, 1, 2)
  75. assert_array_almost_equal(nump.c, [x**2, 2*x + x.conjugate(), 1.0])
  76. assert_array_almost_equal(denomp.c, [x + x.conjugate(), 1.0])
  77. nump, denomp = pade(an, 2, 2)
  78. assert_array_almost_equal(
  79. nump.c,
  80. [x**2 + x*x.conjugate() + x.conjugate()**2, 2*(x + x.conjugate()), 1.0]
  81. )
  82. assert_array_almost_equal(denomp.c, [x.conjugate()**2, x + 2*x.conjugate(), 1.0])