TridiagonalSystemSolver.js 4.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113
  1. import Cartesian3 from "./Cartesian3.js";
  2. import defined from "./defined.js";
  3. import DeveloperError from "./DeveloperError.js";
  4. /**
  5. * Uses the Tridiagonal Matrix Algorithm, also known as the Thomas Algorithm, to solve
  6. * a system of linear equations where the coefficient matrix is a tridiagonal matrix.
  7. *
  8. * @namespace TridiagonalSystemSolver
  9. */
  10. const TridiagonalSystemSolver = {};
  11. /**
  12. * Solves a tridiagonal system of linear equations.
  13. *
  14. * @param {number[]} diagonal An array with length <code>n</code> that contains the diagonal of the coefficient matrix.
  15. * @param {number[]} lower An array with length <code>n - 1</code> that contains the lower diagonal of the coefficient matrix.
  16. * @param {number[]} upper An array with length <code>n - 1</code> that contains the upper diagonal of the coefficient matrix.
  17. * @param {Cartesian3[]} right An array of Cartesians with length <code>n</code> that is the right side of the system of equations.
  18. *
  19. * @exception {DeveloperError} diagonal and right must have the same lengths.
  20. * @exception {DeveloperError} lower and upper must have the same lengths.
  21. * @exception {DeveloperError} lower and upper must be one less than the length of diagonal.
  22. *
  23. * @performance Linear time.
  24. *
  25. * @example
  26. * const lowerDiagonal = [1.0, 1.0, 1.0, 1.0];
  27. * const diagonal = [2.0, 4.0, 4.0, 4.0, 2.0];
  28. * const upperDiagonal = [1.0, 1.0, 1.0, 1.0];
  29. * const rightHandSide = [
  30. * new Cesium.Cartesian3(410757.0, -1595711.0, 1375302.0),
  31. * new Cesium.Cartesian3(-5986705.0, -2190640.0, 1099600.0),
  32. * new Cesium.Cartesian3(-12593180.0, 288588.0, -1755549.0),
  33. * new Cesium.Cartesian3(-5349898.0, 2457005.0, -2685438.0),
  34. * new Cesium.Cartesian3(845820.0, 1573488.0, -1205591.0)
  35. * ];
  36. *
  37. * const solution = Cesium.TridiagonalSystemSolver.solve(lowerDiagonal, diagonal, upperDiagonal, rightHandSide);
  38. *
  39. * @returns {Cartesian3[]} An array of Cartesians with length <code>n</code> that is the solution to the tridiagonal system of equations.
  40. */
  41. TridiagonalSystemSolver.solve = function (lower, diagonal, upper, right) {
  42. //>>includeStart('debug', pragmas.debug);
  43. if (!defined(lower) || !(lower instanceof Array)) {
  44. throw new DeveloperError("The array lower is required.");
  45. }
  46. if (!defined(diagonal) || !(diagonal instanceof Array)) {
  47. throw new DeveloperError("The array diagonal is required.");
  48. }
  49. if (!defined(upper) || !(upper instanceof Array)) {
  50. throw new DeveloperError("The array upper is required.");
  51. }
  52. if (!defined(right) || !(right instanceof Array)) {
  53. throw new DeveloperError("The array right is required.");
  54. }
  55. if (diagonal.length !== right.length) {
  56. throw new DeveloperError("diagonal and right must have the same lengths.");
  57. }
  58. if (lower.length !== upper.length) {
  59. throw new DeveloperError("lower and upper must have the same lengths.");
  60. } else if (lower.length !== diagonal.length - 1) {
  61. throw new DeveloperError(
  62. "lower and upper must be one less than the length of diagonal."
  63. );
  64. }
  65. //>>includeEnd('debug');
  66. const c = new Array(upper.length);
  67. const d = new Array(right.length);
  68. const x = new Array(right.length);
  69. let i;
  70. for (i = 0; i < d.length; i++) {
  71. d[i] = new Cartesian3();
  72. x[i] = new Cartesian3();
  73. }
  74. c[0] = upper[0] / diagonal[0];
  75. d[0] = Cartesian3.multiplyByScalar(right[0], 1.0 / diagonal[0], d[0]);
  76. let scalar;
  77. for (i = 1; i < c.length; ++i) {
  78. scalar = 1.0 / (diagonal[i] - c[i - 1] * lower[i - 1]);
  79. c[i] = upper[i] * scalar;
  80. d[i] = Cartesian3.subtract(
  81. right[i],
  82. Cartesian3.multiplyByScalar(d[i - 1], lower[i - 1], d[i]),
  83. d[i]
  84. );
  85. d[i] = Cartesian3.multiplyByScalar(d[i], scalar, d[i]);
  86. }
  87. scalar = 1.0 / (diagonal[i] - c[i - 1] * lower[i - 1]);
  88. d[i] = Cartesian3.subtract(
  89. right[i],
  90. Cartesian3.multiplyByScalar(d[i - 1], lower[i - 1], d[i]),
  91. d[i]
  92. );
  93. d[i] = Cartesian3.multiplyByScalar(d[i], scalar, d[i]);
  94. x[x.length - 1] = d[d.length - 1];
  95. for (i = x.length - 2; i >= 0; --i) {
  96. x[i] = Cartesian3.subtract(
  97. d[i],
  98. Cartesian3.multiplyByScalar(x[i + 1], c[i], x[i]),
  99. x[i]
  100. );
  101. }
  102. return x;
  103. };
  104. export default TridiagonalSystemSolver;