1""" 2Solution of equations using dense matrices. 3 4The dense matrix is stored as a list of lists. 5 6""" 7 8import copy 9 10from sympy.core.power import isqrt 11from sympy.core.symbol import symbols 12from sympy.matrices.densetools import ( 13 augment, col, conjugate_transpose, eye, rowadd, rowmul) 14from sympy.utilities.exceptions import SymPyDeprecationWarning 15 16SymPyDeprecationWarning( 17 feature="densesolve", 18 issue=12695, 19 deprecated_since_version="1.1").warn() 20 21 22def row_echelon(matlist, K): 23 """ 24 Returns the row echelon form of a matrix with diagonal elements 25 reduced to 1. 26 27 Examples 28 ======== 29 30 >>> from sympy.matrices.densesolve import row_echelon 31 >>> from sympy import QQ 32 >>> a = [ 33 ... [QQ(3), QQ(7), QQ(4)], 34 ... [QQ(2), QQ(4), QQ(5)], 35 ... [QQ(6), QQ(2), QQ(3)]] 36 >>> row_echelon(a, QQ) 37 [[1, 7/3, 4/3], [0, 1, -7/2], [0, 0, 1]] 38 39 See Also 40 ======== 41 42 rref 43 """ 44 result_matlist = copy.deepcopy(matlist) 45 nrow = len(result_matlist) 46 for i in range(nrow): 47 if (result_matlist[i][i] != 1 and result_matlist[i][i] != 0): 48 rowmul(result_matlist, i, 1/result_matlist[i][i], K) 49 for j in range(i + 1, nrow): 50 if (result_matlist[j][i] != 0): 51 rowadd(result_matlist, j, i, -result_matlist[j][i], K) 52 return result_matlist 53 54 55def rref(matlist, K): 56 """ 57 Returns the reduced row echelon form of a Matrix. 58 59 Examples 60 ======== 61 62 >>> from sympy.matrices.densesolve import rref 63 >>> from sympy import QQ 64 >>> a = [ 65 ... [QQ(1), QQ(2), QQ(1)], 66 ... [QQ(-2), QQ(-3), QQ(1)], 67 ... [QQ(3), QQ(5), QQ(0)]] 68 >>> rref(a, QQ) 69 [[1, 0, -5], [0, 1, 3], [0, 0, 0]] 70 71 See Also 72 ======== 73 74 row_echelon 75 """ 76 result_matlist = copy.deepcopy(matlist) 77 result_matlist = row_echelon(result_matlist, K) 78 nrow = len(result_matlist) 79 for i in range(nrow): 80 if result_matlist[i][i] == 1: 81 for j in range(i): 82 rowadd(result_matlist, j, i, -result_matlist[j][i], K) 83 return result_matlist 84 85 86def LU(matlist, K, reverse = 0): 87 """ 88 It computes the LU decomposition of a matrix and returns L and U 89 matrices. 90 91 Examples 92 ======== 93 94 >>> from sympy.matrices.densesolve import LU 95 >>> from sympy import QQ 96 >>> a = [ 97 ... [QQ(1), QQ(2), QQ(3)], 98 ... [QQ(2), QQ(-4), QQ(6)], 99 ... [QQ(3), QQ(-9), QQ(-3)]] 100 >>> LU(a, QQ) 101 ([[1, 0, 0], [2, 1, 0], [3, 15/8, 1]], [[1, 2, 3], [0, -8, 0], [0, 0, -12]]) 102 103 See Also 104 ======== 105 106 upper_triangle 107 lower_triangle 108 """ 109 nrow = len(matlist) 110 new_matlist1, new_matlist2 = eye(nrow, K), copy.deepcopy(matlist) 111 for i in range(nrow): 112 for j in range(i + 1, nrow): 113 if (new_matlist2[j][i] != 0): 114 new_matlist1[j][i] = new_matlist2[j][i]/new_matlist2[i][i] 115 rowadd(new_matlist2, j, i, -new_matlist2[j][i]/new_matlist2[i][i], K) 116 return new_matlist1, new_matlist2 117 118 119def cholesky(matlist, K): 120 """ 121 Performs the cholesky decomposition of a Hermitian matrix and 122 returns L and it's conjugate transpose. 123 124 Examples 125 ======== 126 127 >>> from sympy.matrices.densesolve import cholesky 128 >>> from sympy import QQ 129 >>> cholesky([[QQ(25), QQ(15), QQ(-5)], [QQ(15), QQ(18), QQ(0)], [QQ(-5), QQ(0), QQ(11)]], QQ) 130 ([[5, 0, 0], [3, 3, 0], [-1, 1, 3]], [[5, 3, -1], [0, 3, 1], [0, 0, 3]]) 131 132 See Also 133 ======== 134 135 cholesky_solve 136 """ 137 new_matlist = copy.deepcopy(matlist) 138 nrow = len(new_matlist) 139 L = eye(nrow, K) 140 for i in range(nrow): 141 for j in range(i + 1): 142 a = K.zero 143 for k in range(j): 144 a += L[i][k]*L[j][k] 145 if i == j: 146 L[i][j] = isqrt(new_matlist[i][j] - a) 147 else: 148 L[i][j] = (new_matlist[i][j] - a)/L[j][j] 149 return L, conjugate_transpose(L, K) 150 151 152def LDL(matlist, K): 153 """ 154 Performs the LDL decomposition of a hermitian matrix and returns L, D and 155 transpose of L. Only applicable to rational entries. 156 157 Examples 158 ======== 159 160 >>> from sympy.matrices.densesolve import LDL 161 >>> from sympy import QQ 162 163 >>> a = [ 164 ... [QQ(4), QQ(12), QQ(-16)], 165 ... [QQ(12), QQ(37), QQ(-43)], 166 ... [QQ(-16), QQ(-43), QQ(98)]] 167 >>> LDL(a, QQ) 168 ([[1, 0, 0], [3, 1, 0], [-4, 5, 1]], [[4, 0, 0], [0, 1, 0], [0, 0, 9]], [[1, 3, -4], [0, 1, 5], [0, 0, 1]]) 169 170 """ 171 new_matlist = copy.deepcopy(matlist) 172 nrow = len(new_matlist) 173 L, D = eye(nrow, K), eye(nrow, K) 174 for i in range(nrow): 175 for j in range(i + 1): 176 a = K.zero 177 for k in range(j): 178 a += L[i][k]*L[j][k]*D[k][k] 179 if i == j: 180 D[j][j] = new_matlist[j][j] - a 181 else: 182 L[i][j] = (new_matlist[i][j] - a)/D[j][j] 183 return L, D, conjugate_transpose(L, K) 184 185 186def upper_triangle(matlist, K): 187 """ 188 Transforms a given matrix to an upper triangle matrix by performing 189 row operations on it. 190 191 Examples 192 ======== 193 194 >>> from sympy.matrices.densesolve import upper_triangle 195 >>> from sympy import QQ 196 >>> a = [ 197 ... [QQ(4,1), QQ(12,1), QQ(-16,1)], 198 ... [QQ(12,1), QQ(37,1), QQ(-43,1)], 199 ... [QQ(-16,1), QQ(-43,1), QQ(98,1)]] 200 >>> upper_triangle(a, QQ) 201 [[4, 12, -16], [0, 1, 5], [0, 0, 9]] 202 203 See Also 204 ======== 205 206 LU 207 """ 208 copy_matlist = copy.deepcopy(matlist) 209 lower_triangle, upper_triangle = LU(copy_matlist, K) 210 return upper_triangle 211 212 213def lower_triangle(matlist, K): 214 """ 215 Transforms a given matrix to a lower triangle matrix by performing 216 row operations on it. 217 218 Examples 219 ======== 220 221 >>> from sympy.matrices.densesolve import lower_triangle 222 >>> from sympy import QQ 223 >>> a = [ 224 ... [QQ(4,1), QQ(12,1), QQ(-16)], 225 ... [QQ(12,1), QQ(37,1), QQ(-43,1)], 226 ... [QQ(-16,1), QQ(-43,1), QQ(98,1)]] 227 >>> lower_triangle(a, QQ) 228 [[1, 0, 0], [3, 1, 0], [-4, 5, 1]] 229 230 See Also 231 ======== 232 233 LU 234 """ 235 copy_matlist = copy.deepcopy(matlist) 236 lower_triangle, upper_triangle = LU(copy_matlist, K, reverse = 1) 237 return lower_triangle 238 239 240def rref_solve(matlist, variable, constant, K): 241 """ 242 Solves a system of equations using reduced row echelon form given 243 a matrix of coefficients, a vector of variables and a vector of constants. 244 245 Examples 246 ======== 247 248 >>> from sympy.matrices.densesolve import rref_solve 249 >>> from sympy import QQ 250 >>> from sympy import Dummy 251 >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z') 252 >>> coefficients = [ 253 ... [QQ(25), QQ(15), QQ(-5)], 254 ... [QQ(15), QQ(18), QQ(0)], 255 ... [QQ(-5), QQ(0), QQ(11)]] 256 >>> constants = [ 257 ... [QQ(2)], 258 ... [QQ(3)], 259 ... [QQ(1)]] 260 >>> variables = [ 261 ... [x], 262 ... [y], 263 ... [z]] 264 >>> rref_solve(coefficients, variables, constants, QQ) 265 [[-1/225], [23/135], [4/45]] 266 267 See Also 268 ======== 269 270 row_echelon 271 augment 272 """ 273 new_matlist = copy.deepcopy(matlist) 274 augmented = augment(new_matlist, constant, K) 275 solution = rref(augmented, K) 276 return col(solution, -1) 277 278 279def LU_solve(matlist, variable, constant, K): 280 """ 281 Solves a system of equations using LU decomposition given a matrix 282 of coefficients, a vector of variables and a vector of constants. 283 284 Examples 285 ======== 286 287 >>> from sympy.matrices.densesolve import LU_solve 288 >>> from sympy import QQ 289 >>> from sympy import Dummy 290 >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z') 291 >>> coefficients = [ 292 ... [QQ(2), QQ(-1), QQ(-2)], 293 ... [QQ(-4), QQ(6), QQ(3)], 294 ... [QQ(-4), QQ(-2), QQ(8)]] 295 >>> variables = [ 296 ... [x], 297 ... [y], 298 ... [z]] 299 >>> constants = [ 300 ... [QQ(-1)], 301 ... [QQ(13)], 302 ... [QQ(-6)]] 303 >>> LU_solve(coefficients, variables, constants, QQ) 304 [[2], [3], [1]] 305 306 See Also 307 ======== 308 309 LU 310 forward_substitution 311 backward_substitution 312 """ 313 new_matlist = copy.deepcopy(matlist) 314 nrow = len(new_matlist) 315 L, U = LU(new_matlist, K) 316 y = [[i] for i in symbols('y:%i' % nrow)] 317 forward_substitution(L, y, constant, K) 318 backward_substitution(U, variable, y, K) 319 return variable 320 321 322def cholesky_solve(matlist, variable, constant, K): 323 """ 324 Solves a system of equations using Cholesky decomposition given 325 a matrix of coefficients, a vector of variables and a vector of constants. 326 327 Examples 328 ======== 329 330 >>> from sympy.matrices.densesolve import cholesky_solve 331 >>> from sympy import QQ 332 >>> from sympy import Dummy 333 >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z') 334 >>> coefficients = [ 335 ... [QQ(25), QQ(15), QQ(-5)], 336 ... [QQ(15), QQ(18), QQ(0)], 337 ... [QQ(-5), QQ(0), QQ(11)]] 338 >>> variables = [ 339 ... [x], 340 ... [y], 341 ... [z]] 342 >>> coefficients = [ 343 ... [QQ(2)], 344 ... [QQ(3)], 345 ... [QQ(1)]] 346 >>> cholesky_solve([[QQ(25), QQ(15), QQ(-5)], [QQ(15), QQ(18), QQ(0)], [QQ(-5), QQ(0), QQ(11)]], [[x], [y], [z]], [[QQ(2)], [QQ(3)], [QQ(1)]], QQ) 347 [[-1/225], [23/135], [4/45]] 348 349 See Also 350 ======== 351 352 cholesky 353 forward_substitution 354 backward_substitution 355 """ 356 new_matlist = copy.deepcopy(matlist) 357 nrow = len(new_matlist) 358 L, Lstar = cholesky(new_matlist, K) 359 y = [[i] for i in symbols('y:%i' % nrow)] 360 forward_substitution(L, y, constant, K) 361 backward_substitution(Lstar, variable, y, K) 362 return variable 363 364 365def forward_substitution(lower_triangle, variable, constant, K): 366 """ 367 Performs forward substitution given a lower triangular matrix, a 368 vector of variables and a vector of constants. 369 370 Examples 371 ======== 372 373 >>> from sympy.matrices.densesolve import forward_substitution 374 >>> from sympy import QQ 375 >>> from sympy import Dummy 376 >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z') 377 >>> a = [ 378 ... [QQ(1), QQ(0), QQ(0)], 379 ... [QQ(-2), QQ(1), QQ(0)], 380 ... [QQ(-2), QQ(-1), QQ(1)]] 381 >>> variables = [ 382 ... [x], 383 ... [y], 384 ... [z]] 385 >>> constants = [ 386 ... [QQ(-1)], 387 ... [QQ(13)], 388 ... [QQ(-6)]] 389 >>> forward_substitution(a, variables, constants, QQ) 390 [[-1], [11], [3]] 391 392 See Also 393 ======== 394 395 LU_solve 396 cholesky_solve 397 """ 398 copy_lower_triangle = copy.deepcopy(lower_triangle) 399 nrow = len(copy_lower_triangle) 400 for i in range(nrow): 401 a = K.zero 402 for j in range(i): 403 a += copy_lower_triangle[i][j]*variable[j][0] 404 variable[i][0] = (constant[i][0] - a)/copy_lower_triangle[i][i] 405 return variable 406 407 408def backward_substitution(upper_triangle, variable, constant, K): 409 """ 410 Performs forward substitution given a lower triangular matrix, 411 a vector of variables and a vector constants. 412 413 Examples 414 ======== 415 416 >>> from sympy.matrices.densesolve import backward_substitution 417 >>> from sympy import QQ 418 >>> from sympy import Dummy 419 >>> x, y, z = Dummy('x'), Dummy('y'), Dummy('z') 420 >>> a = [ 421 ... [QQ(2), QQ(-1), QQ(-2)], 422 ... [QQ(0), QQ(4), QQ(-1)], 423 ... [QQ(0), QQ(0), QQ(3)]] 424 >>> variables = [ 425 ... [x], 426 ... [y], 427 ... [z]] 428 >>> constants = [ 429 ... [QQ(-1)], 430 ... [QQ(11)], 431 ... [QQ(3)]] 432 >>> backward_substitution(a, variables, constants, QQ) 433 [[2], [3], [1]] 434 435 See Also 436 ======== 437 438 LU_solve 439 cholesky_solve 440 """ 441 copy_upper_triangle = copy.deepcopy(upper_triangle) 442 nrow = len(copy_upper_triangle) 443 for i in reversed(range(nrow)): 444 a = K.zero 445 for j in reversed(range(i + 1, nrow)): 446 a += copy_upper_triangle[i][j]*variable[j][0] 447 variable[i][0] = (constant[i][0] - a)/copy_upper_triangle[i][i] 448 return variable 449