Cramer's rule: Difference between revisions

Content added Content deleted

Inline

Revision as of 09:47, 12 April 2016

In linear algebra, Cramer's rule is an explicit formula for the solution of a system of linear equations with as many equations as unknowns, valid whenever the system has a unique solution. It expresses the solution in terms of the determinants of the (square) coefficient matrix and of matrices obtained from it by replacing one column by the vector of right hand sides of the equations.

Given

\left\{{\begin{matrix}a_{1}x+b_{1}y+c_{1}z&={\color {red}d_{1}}\\a_{2}x+b_{2}y+c_{2}z&={\color {red}d_{2}}\\a_{3}x+b_{3}y+c_{3}z&={\color {red}d_{3}}\end{matrix}}\right.

which in matrix format is

${\begin{bmatrix}a_{1}&b_{1}&c_{1}\\a_{2}&b_{2}&c_{2}\\a_{3}&b_{3}&c_{3}\end{bmatrix}}{\begin{bmatrix}x\\y\\z\end{bmatrix}}={\begin{bmatrix}{\color {red}d_{1}}\\{\color {red}d_{2}}\\{\color {red}d_{3}}\end{bmatrix}}.$

Then the values of $x,y$ and $z$ can be found as follows:

x={\frac {\begin{vmatrix}{\color {red}d_{1}}&b_{1}&c_{1}\\{\color {red}d_{2}}&b_{2}&c_{2}\\{\color {red}d_{3}}&b_{3}&c_{3}\end{vmatrix}}{\begin{vmatrix}a_{1}&b_{1}&c_{1}\\a_{2}&b_{2}&c_{2}\\a_{3}&b_{3}&c_{3}\end{vmatrix}}},\quad y={\frac {\begin{vmatrix}a_{1}&{\color {red}d_{1}}&c_{1}\\a_{2}&{\color {red}d_{2}}&c_{2}\\a_{3}&{\color {red}d_{3}}&c_{3}\end{vmatrix}}{\begin{vmatrix}a_{1}&b_{1}&c_{1}\\a_{2}&b_{2}&c_{2}\\a_{3}&b_{3}&c_{3}\end{vmatrix}}},{\text{ and }}z={\frac {\begin{vmatrix}a_{1}&b_{1}&{\color {red}d_{1}}\\a_{2}&b_{2}&{\color {red}d_{2}}\\a_{3}&b_{3}&{\color {red}d_{3}}\end{vmatrix}}{\begin{vmatrix}a_{1}&b_{1}&c_{1}\\a_{2}&b_{2}&c_{2}\\a_{3}&b_{3}&c_{3}\end{vmatrix}}}.

Task

Given the following system of equations:

{\begin{cases}2w-x+5y+z=-3\\3w+2x+2y-6z=-32\\w+3x+3y-z=-47\\5w-2x-3y+3z=49\\\end{cases}}

solve for $w$ , $x$ , $y$ and $z$ , using Cramer's rule.

EchoLisp

<lang scheme> (lib 'matrix) (string-delimiter "") (define (cramer A B (X)) ;; --> vector X (define ∆ (determinant A)) (for/vector [(i (matrix-col-num A))] (set! X (matrix-set-col! (array-copy A) i B)) (// (determinant X) ∆)))

(define (task) (define A (list->array

 	'( 2 -1 5 1 3 2 2 -6 1 3 3 -1 5 -2 -3 3) 4 4))

(define B #(-3 -32 -47 49)) (writeln "Solving A * X = B") (array-print A) (writeln "B = " B) (writeln "X = " (cramer A B))) </lang>

Output:

(task)
Solving A * X = B    
  2   -1   5    1  
  3   2    2    -6 
  1   3    3    -1 
  5   -2   -3   3  
B =      #( -3 -32 -47 49)    
X =      #( 2 -12 -4 1)

Fortran

In Numerical Methods That Work (Usually), in the section What not to compute, F. S. Acton remarks "...perhaps we should be glad he didn't resort to Cramer's rule (still taught as the practical method in some high schools) and solve his equations as the ratios of determinants - a process that requires labor proportional to N! if done in the schoolboy manner. The contrast with N³ can be startling! Further on, "Having hinted darkly at my computational fundamentalism, it is probably time to commit to a public heresy by denouncing recursive calculations. I have never seen a numerical problem arising from the physical world that was best calculated by a recursive subroutine..."

Since this problem requires use of Cramer's rule, one might as well be hung for a sheep instead of a lamb, so the traditions of Old Fortran will be ignored and the fearsome RECURSIVE specification employed so that the determinants will be calculated recursively, all the way down to N = 1 even though the N = 2 case is easy. This requires F90 and later.

     MODULE BAD IDEA	!Employ Cramer's rule for solving A.x = b...
      INTEGER MSG	!Might as well be here.
      CONTAINS		!The details.
       SUBROUTINE SHOWMATRIX(A)	!With nice vertical bars.
        DOUBLE PRECISION A(:,:)	!The matrix.
        INTEGER R,N			!Assistants.
         N = SIZE(A, DIM = 1)		!Instead of passing an explicit parameter.
         DO R = 1,N			!Work down the rows.
           WRITE (MSG,1) A(R,:)		!Fling forth a row at a go.
   1       FORMAT ("|",<N>F12.3,"|")		!Bounded by bars.
         END DO			!On to the next row.
       END SUBROUTINE SHOWMATRIX	!Furrytran's default order is the transpose.

       RECURSIVE DOUBLE PRECISION FUNCTION DET(A)	!Determine the determinant.
        DOUBLE PRECISION A(:,:)	!The square matrix, order N.
        DOUBLE PRECISION MINOR(SIZE(A,DIM=1) - 1,SIZE(A,DIM=1) - 1)	!Order N - 1.
        DOUBLE PRECISION D		!A waystation.
        INTEGER C,N			!Steppers.
         N = SIZE(A, DIM = 1)		!Suplied via secret parameters.
         IF (N .LE. 0) THEN		!This really ought not happen.
           STOP "DET: null array!"		!But I'm endlessly suspicious.
         ELSE IF (N .NE. SIZE(A, DIM = 2)) THEN	!And I'd rather have a decent message
           STOP "DET: not a square array!"	!In place of a crashed run.
         ELSE IF (N .EQ. 1) THEN	!Alright, now get on with it.
           DET = A(1,1)		!This is really easy.
         ELSE				!But otherwise,
           D = 0			!Here we go.
           DO C = 1,N			!Step along the columns of the first row.
             CALL FILLMINOR(C)			!Produce the auxiliary array for each column.
             IF (MOD(C,2) .EQ. 0) THEN		!Odd or even case?
               D = D - A(1,C)*DET(MINOR)	!Even: subtract.
              ELSE				!Otherwise,
               D = D + A(1,C)*DET(MINOR)	!Odd: add.
             END IF				!So much for that term.
           END DO			!On to the next.
           DET = D		!Declare the result.
         END IF	!So much for that.
        CONTAINS	!An assistant.
          SUBROUTINE FILLMINOR(CC)	!Corresponding to A(1,CC).
          INTEGER CC	!The column being omitted.
          INTEGER R	!A stepper.
          DO R = 2,N		!Ignoring the first row,
            MINOR(R - 1,1:CC - 1) = A(R,1:CC - 1)	!Copy columns 1 to CC - 1. There may be none.
            MINOR(R - 1,CC:) = A(R,CC + 1:)		!And from CC + 1 to N. There may be none.
          END DO		!On to the next row.
         END SUBROUTINE FILLMINOR	!Divide and conquer.
       END FUNCTION DET	!Rather than mess with permutations.

       SUBROUTINE CRAMER(A,X,B)	!Solve A.x= b, where A is a matrix...

Careful! The array must be A(N,N), and not say A(100,100) of which only up to N = 6 are in use.

        DOUBLE PRECISION A(:,:)	!A square matrix. I hope.
        DOUBLE PRECISION B(:),X(:)	!Column matrices look rather like 1-D arrays.
        DOUBLE PRECISION AUX(SIZE(A,DIM=1),SIZE(A,DIM=1))	!Can't say "LIKE A", as in pl/i, alas.
        DOUBLE PRECISION D	!To be calculated once.
        INTEGER C		!A stepper.
         N = SIZE(A, DIM = 1)	!Alright, what's the order of battle?
         D = DET(A)		!Once only.
         IF (D.EQ.0) STOP "Cramer: zero determinant!"	!Surely, this won't happen...
         AUX = A		!Prepare the assistant.
         DO C = 1,N		!Step across the columns.
           IF (C.GT.1) AUX(1:N,C - 1) = A(1:N,C - 1)	!Repair previous damage.
           AUX(1:N,C) = B(1:N)		!Place the current damage.
           X(C) = DET(AUX)/D		!The result!
         END DO		!On to the next column.
       END SUBROUTINE CRAMER	!This looks really easy!
     END MODULE BAD IDEA	!But actually, it is a bad idea for N > 2.

     PROGRAM TEST	!Try it and see.
     USE BAD IDEA	!Just so.
     DOUBLE PRECISION, ALLOCATABLE ::A(:,:), B(:), X(:)	!Custom work areas.
     INTEGER N,R	!Assistants..
     INTEGER INF	!An I/O unit.

     MSG = 6		!Output.
     INF = 10		!For reading test data.
     OPEN (INF,FILE="Test.dat",STATUS="OLD",ACTION="READ")	!As in this file..

Chew into the next problem.

  10 IF (ALLOCATED(A)) DEALLOCATE(A)	!First,
     IF (ALLOCATED(B)) DEALLOCATE(B)	!Get rid of
     IF (ALLOCATED(X)) DEALLOCATE(X)	!The hired help.
     READ (INF,*,END = 100) N		!Since there is a new order.
     IF (N.LE.0) GO TO 100		!Perhaps a final order.
     WRITE (MSG,11) N			!Othewise, announce prior to acting.
  11 FORMAT ("Order ",I0," matrix as follows...")	!In case something goes wrong.
     ALLOCATE(A(N,N))			!For instance,
     ALLOCATE(B(N))			!Out of memory.
     ALLOCATE(X(N))			!But otherwise, a tailored fit.
     DO R = 1,N			!Now read in the values for the matrix.
       READ(INF,*, END = 666) A(R,:)	!One row at a go.
     END DO				!In free format.
     CALL SHOWMATRIX(A)		!Show what we have managed to obtain.
     READ (INF,*) B			!Now get the b of A.x = b.
     WRITE (MSG,12) B			!Reveal.
  12 FORMAT ("In the scheme A.x = b, b = ",<N>F12.6)	!In case something goes wrong.
     CALL CRAMER(A,X,B)		!The deed!
     WRITE (MSG,13) X			!The result!
  13 FORMAT ("    Via Cramer's rule, x = ",<N>F12.6)	!Aligned with FORMAT 12.
     GO TO 10				!And try for another test problem.

Complaints.

 666 WRITE (MSG,667) R		!I know where I came from.
 667 FORMAT ("End-of-file while reading row ",I0,"!")	!So, explain.

Closedown.

 100 WRITE (6,*) "That was interesting."	!Quite.
     END	!Open files are closed, allocated memory is released.

</lang> Output: file Test.dat has been filled with the numbers of the example.

Order 4 matrix as follows...
|       2.000      -1.000       5.000       1.000|
|       3.000       2.000       2.000      -6.000|
|       1.000       3.000       3.000      -1.000|
|       5.000      -2.000      -3.000       3.000|
In the scheme A.x = b, b =    -3.000000  -32.000000  -47.000000   49.000000
    Via Cramer's rule, x =     2.000000  -12.000000   -4.000000    1.000000
 That was interesting.

J

Implementation:

<lang J>cramer=:4 :0

 A=. x [ b=. y
 det=. -/ .*
 A %~&det (i.#A) b"_`[`]}&.|:"0 2 A

)</lang>

Task data:

<lang J>A=: _&".;._2]t=: 0 :0

)

b=: _3 _32 _47 49</lang>

Task example:

<lang J> A cramer b 2 _12 _4 1</lang>

Perl

<lang perl>use Math::Matrix;

sub cramers_rule {

   my ($A, $terms) = @_;
   my @solutions;
   my $det = $A->determinant;
   foreach my $i (0 .. $#{$A}) {
       my $Ai = $A->clone;
       foreach my $j (0 .. $#{$terms}) {
           $Ai->[$j][$i] = $terms->[$j];
       }
       push @solutions, $Ai->determinant / $det;
   }
   @solutions;

}

my $matrix = Math::Matrix->new(

   [2, -1,  5,  1],
   [3,  2,  2, -6],
   [1,  3,  3, -1],
   [5, -2, -3,  3],

);

my $free_terms = [-3, -32, -47, 49]; my ($w, $x, $y, $z) = cramers_rule($matrix, $free_terms);

print "w = $w\n"; print "x = $x\n"; print "y = $y\n"; print "z = $z\n";</lang>

Output:

w = 2
x = -12
y = -4
z = 1

Perl 6

<lang perl6>sub det(@matrix) {

   my @a = @matrix.map: { [|$_] };
   my $sign = +1;
   my $pivot = 1;
   for ^@a -> $k {
     my @r = ($k+1 .. @a.end);
     my $previous-pivot = $pivot;
     if 0 == ($pivot = @a[$k][$k]) {
       (my $s = @r.first: { @a[$_][$k] != 0 }) // return 0;
       (@a[$s],@a[$k]) = (@a[$k], @a[$s]);
       my $pivot = @a[$k][$k];
       $sign = -$sign;
     }
     for @r X @r -> ($i, $j) {
       ((@a[$i][$j] *= $pivot) -= @a[$i][$k]*@a[$k][$j]) /= $previous-pivot;
     }
   }
   $sign * $pivot

}

sub cramers_rule(@A, @terms) {

   gather for ^@A -> $i {
       my @Ai = @A.map: { [|$_] };
       for ^@terms -> $j {
           @Ai[$j][$i] = @terms[$j];
       }
       take det(@Ai);
   } »/» det(@A);

}

my @matrix = (

   [2, -1,  5,  1],
   [3,  2,  2, -6],
   [1,  3,  3, -1],
   [5, -2, -3,  3],

);

my @free_terms = (-3, -32, -47, 49); my ($w, $x, $y, $z) = |cramers_rule(@matrix, @free_terms);

say "w = $w"; say "x = $x"; say "y = $y"; say "z = $z";</lang>

Output:

w = 2
x = -12
y = -4
z = 1

Ruby

<lang ruby>require 'matrix'

def cramers_rule(a, terms)

   solutions = []
   det = a.determinant
   size = a.row_size;
   size.times do |i|
       cols = []
       (0..i-1).each { |j|
           cols << a.column(j)
       }
       cols << terms
       (i+1...size).each { |j|
           cols << a.column(j)
       }
       ai = Matrix.columns(cols)
       solutions << ai.determinant / det
   end
   solutions

end

matrix = Matrix[

   [2, -1,  5,  1],
   [3,  2,  2, -6],
   [1,  3,  3, -1],
   [5, -2, -3,  3],

]

vector = [-3, -32, -47, 49] puts cramers_rule(matrix, vector)</lang>

Output:

Sidef

<lang ruby>func det(a) {

   a = a.map{.map{_}}
   var sign = +1
   var pivot = 1

   for k in ^a {
     var r = (k+1 .. a.end)
     var previous_pivot = pivot

     if ((pivot = a[k][k])==0) {
       a.swap(r.first_by { |i|
           a[i][k] != 0
       } \\ (return 0), k)
       pivot = a[k][k]
       sign = -sign
     }

     for i,j in (r ~X r) {
       ((a[i][j] *= pivot)
           -= a[i][k]*a[k][j]
       ) /= previous_pivot
     }
   }
   sign * pivot

}

func cramers_rule(A, terms) {

   gather {
       for i in ^A {
           var Ai = A.map{.map{_}}
           for j in ^terms {
               Ai[j][i] = terms[j]
           }
           take(det(Ai))
       }
   } »/» det(A)

}

var matrix = [

   [2, -1,  5,  1],
   [3,  2,  2, -6],
   [1,  3,  3, -1],
   [5, -2, -3,  3],

]

var free_terms = [-3, -32, -47, 49] var (w, x, y, z) = cramers_rule(matrix, free_terms)...

say "w = #{w}" say "x = #{x}" say "y = #{y}" say "z = #{z}"</lang>

Output:

w = 2
x = -12
y = -4
z = 1