next | previous | forward | backward | up | top | index | toc | Macaulay2 web site
Macaulay2Doc :: solve

solve -- solve linear equation(s)

Synopsis

Description

i1 : kk = ZZ/101;
i2 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk

o2 = | 1  2 3  4  |
     | 1  3 6  10 |
     | 19 7 11 13 |

              3        4
o2 : Matrix kk  <--- kk
i3 : b = matrix"1;1;1" ** kk

o3 = | 1 |
     | 1 |
     | 1 |

              3        1
o3 : Matrix kk  <--- kk
i4 : x = solve(A,b)

o4 = | 2  |
     | -1 |
     | 34 |
     | 0  |

              4        1
o4 : Matrix kk  <--- kk
i5 : A*x-b

o5 = 0

              3        1
o5 : Matrix kk  <--- kk
i6 : kk = GF(25)

o6 = kk

o6 : GaloisField
i7 : a = kk_0

o7 = a

o7 : kk
i8 : A = matrix"a,a+1,a+2,3a,4;a-1,1,2a,6,10;19,7,a,11,13" ** kk

o8 = | a   a+1 a+2 -2a -1 |
     | a-1 1   2a  1   0  |
     | -1  2   a   1   -2 |

              3        5
o8 : Matrix kk  <--- kk
i9 : b = matrix"1;-a+1;1" ** kk

o9 = | 1    |
     | -a+1 |
     | 1    |

              3        1
o9 : Matrix kk  <--- kk
i10 : x = solve(A,b)

o10 = | -a    |
      | -2a+1 |
      | -2a   |
      | 0     |
      | 0     |

               5        1
o10 : Matrix kk  <--- kk
i11 : A*x-b

o11 = 0

               3        1
o11 : Matrix kk  <--- kk
i12 : kk = QQ

o12 = QQ

o12 : Ring
i13 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk

o13 = | 1  2 3  4  |
      | 1  3 6  10 |
      | 19 7 11 13 |

               3        4
o13 : Matrix QQ  <--- QQ
i14 : b = matrix"1;1;1" ** kk

o14 = | 1 |
      | 1 |
      | 1 |

               3        1
o14 : Matrix QQ  <--- QQ
i15 : x = solve(A,b)

o15 = | -7/47  |
      | 54/47  |
      | -18/47 |
      | 0      |

               4        1
o15 : Matrix QQ  <--- QQ
i16 : A*x-b

o16 = 0

               3        1
o16 : Matrix QQ  <--- QQ

Over RR53 or CC53, if the matrix A is non-singular and square, then highly optimized lapack routines will be called.

i17 : printingPrecision = 4;
i18 : A = matrix "1,2,3;1,3,6;19,7,11" ** RR

o18 = | 1  2 3  |
      | 1  3 6  |
      | 19 7 11 |

                 3          3
o18 : Matrix RR    <--- RR
               53         53
i19 : b = matrix "1;1;1" ** RR

o19 = | 1 |
      | 1 |
      | 1 |

                 3          1
o19 : Matrix RR    <--- RR
               53         53
i20 : x = solve(A,b)

o20 = | -.1489 |
      | 1.149  |
      | -.383  |

                 3          1
o20 : Matrix RR    <--- RR
               53         53
i21 : A*x-b

o21 = | 2.22e-16  |
      | -2.22e-16 |
      | 0         |

                 3          1
o21 : Matrix RR    <--- RR
               53         53
i22 : norm oo

o22 = 2.22044604925031e-16

o22 : RR (of precision 53)
i23 : clean(1e-15, A*x-b)

o23 = | 0 |
      | 0 |
      | 0 |

                 3          1
o23 : Matrix RR    <--- RR
               53         53

If you know that your matrix is square, and invertible, then providing the hint: MaximalRank=>true allows Macaulay2 to choose the fastest routines. For small matrix sizes, it should not be too noticeable, but for large matrices, the difference in time taken can be dramatic.

i24 : printingPrecision = 4;
i25 : N = 40

o25 = 40
i26 : A = mutableMatrix(CC_53, N, N); fillMatrix A;
i28 : B = mutableMatrix(CC_53, N, 2); fillMatrix B;
i30 : time X = solve(A,B);
     -- used 0.000280063 seconds
i31 : time X = solve(A,B, MaximalRank=>true);
     -- used 0.000113046 seconds
i32 : norm(A*X-B)

o32 = 5.11185069084045e-15

o32 : RR (of precision 53)

Over higher precision RR or CC, these routines will be much slower than the lower precision lapack routines.

i33 : N = 100

o33 = 100
i34 : A = mutableMatrix(CC_100, N, N); fillMatrix A;
i36 : B = mutableMatrix(CC_100, N, 2); fillMatrix B;
i38 : time X = solve(A,B);
     -- used 0.189148 seconds
i39 : time X = solve(A,B, MaximalRank=>true);
     -- used 0.188485 seconds
i40 : norm(A*X-B)

o40 = 1.49157827468970981408235588593e-28

o40 : RR (of precision 100)

Giving the option ClosestFit=>true, in the case when the field is RR or CC, uses a least squares algorithm to find a best fit solution.

i41 : kk = RR_53;
i42 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk

o42 = | 1  2 3  4  |
      | 1  3 6  10 |
      | 19 7 11 13 |

                 3          4
o42 : Matrix RR    <--- RR
               53         53
i43 : b = matrix"1;1;1" ** kk

o43 = | 1 |
      | 1 |
      | 1 |

                 3          1
o43 : Matrix RR    <--- RR
               53         53
i44 : x1 = solve(A,b, ClosestFit=>true)

o44 = | -.1899 |
      | .6399  |
      | .3367  |
      | -.275  |

                 4          1
o44 : Matrix RR    <--- RR
               53         53
i45 : A*x1-b

o45 = | -5.551e-16 |
      | -3.331e-16 |
      | 0          |

                 3          1
o45 : Matrix RR    <--- RR
               53         53

Giving both options ClosestFit and MaximalRank allows Macaulay2 to call a faster algorithm.

i46 : x2 = solve(A,b, ClosestFit=>true, MaximalRank=>true)

o46 = | -.1899 |
      | .6399  |
      | .3367  |
      | -.275  |

                 4          1
o46 : Matrix RR    <--- RR
               53         53
i47 : A*x2-b

o47 = | 0         |
      | 0         |
      | 3.109e-15 |

                 3          1
o47 : Matrix RR    <--- RR
               53         53

Caveat

(1) This function is limited in scope, but has been designed to be much faster than generic algorithms. (2) If the matrix is a square invertible matrix, giving the option MaximalRank=>true can strongly speed up the computation. (3) For mutable matrices, this function is only currently implemented for densely encoded matrices.

See also

Ways to use solve :