linear algebra
Cramer's rule
system of linear equations


Given


 




{





a

1


x
+

b

1


y
+

c

1


z


=



d

1









a

2


x
+

b

2


y
+

c

2


z


=



d

2









a

3


x
+

b

3


y
+

c

3


z


=



d

3












{\displaystyle \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



 





[




a

1





b

1





c

1







a

2





b

2





c

2







a

3





b

3





c

3





]




[



x




y




z



]


=


[






d

1











d

2











d

3







]


.


{\displaystyle {\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


{\displaystyle x,y}

 and 



z


{\displaystyle z}

 can be found as follows:







x
=



|






d

1







b

1





c

1









d

2







b

2





c

2









d

3







b

3





c

3





|


|




a

1





b

1





c

1







a

2





b

2





c

2







a

3





b

3





c

3





|



,

y
=



|




a

1







d

1







c

1







a

2







d

2







c

2







a

3







d

3







c

3





|


|




a

1





b

1





c

1







a

2





b

2





c

2







a

3





b

3





c

3





|



,

 and 

z
=



|




a

1





b

1







d

1









a

2





b

2







d

2









a

3





b

3







d

3







|


|




a

1





b

1





c

1







a

2





b

2





c

2







a

3





b

3





c

3





|



.


{\displaystyle 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:









{



2
w
−
x
+
5
y
+
z
=
−
3




3
w
+
2
x
+
2
y
−
6
z
=
−
32




w
+
3
x
+
3
y
−
z
=
−
47




5
w
−
2
x
−
3
y
+
3
z
=
49








{\displaystyle {\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


{\displaystyle w}

, 



x


{\displaystyle x}

, 



y


{\displaystyle y}

 and 



z


{\displaystyle z}

, using Cramer's rule.




use std::ops::{Index, IndexMut};
 
fn main() {
    let m = matrix(
        vec![
            2., -1., 5., 1., 3., 2., 2., -6., 1., 3., 3., -1., 5., -2., -3., 3.,
        ],
        4,
    );
    let mm = m.solve(&vec![-3., -32., -47., 49.]);
    println!("{:?}", mm);
}
 
#[derive(Clone)]
struct Matrix {
    elts: Vec<f64>,
    dim: usize,
}
 
impl Matrix {
    // Compute determinant using cofactor method
    // Using Gaussian elimination would have been more efficient, but it also solves the linear
    // system, so…
    fn det(&self) -> f64 {
        match self.dim {
            0 => 0.,
            1 => self[0][0],
            2 => self[0][0] * self[1][1] - self[0][1] * self[1][0],
            d => {
                let mut acc = 0.;
                let mut signature = 1.;
                for k in 0..d {
                    acc += signature * self[0][k] * self.comatrix(0, k).det();
                    signature *= -1.
                }
                acc
            }
        }
    }
 
    // Solve linear systems using Cramer's method
    fn solve(&self, target: &Vec<f64>) -> Vec<f64> {
        let mut solution: Vec<f64> = vec![0.; self.dim];
        let denominator = self.det();
        for j in 0..self.dim {
            let mut mm = self.clone();
            for i in 0..self.dim {
                mm[i][j] = target[i]
            }
            solution[j] = mm.det() / denominator
        }
        solution
    }
 
    // Compute the cofactor matrix for determinant computations
    fn comatrix(&self, k: usize, l: usize) -> Matrix {
        let mut v: Vec<f64> = vec![];
        for i in 0..self.dim {
            for j in 0..self.dim {
                if i != k && j != l {
                    v.push(self[i][j])
                }
            }
        }
        matrix(v, self.dim - 1)
    }
}
 
fn matrix(elts: Vec<f64>, dim: usize) -> Matrix {
    assert_eq!(elts.len(), dim * dim);
    Matrix { elts, dim }
}
 
impl Index<usize> for Matrix {
    type Output = [f64];
 
    fn index(&self, i: usize) -> &Self::Output {
        let m = self.dim;
        &self.elts[m * i..m * (i + 1)]
    }
}
 
impl IndexMut<usize> for Matrix {
    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
        let m = self.dim;
        &mut self.elts[m * i..m * (i + 1)]
    }
}