linear_systems/
linear_systems.rs

1// cargo run --example linear_systems
2
3extern crate ohsl;
4
5pub use ohsl::vector::{Vector, Vec64};
6pub use ohsl::matrix::{Matrix, Mat64};
7pub use ohsl::complex::{Complex, Cmplx};
8
9
10
11fn main() {
12    println!( "------------------ Linear system ------------------" );
13    println!( "  * Solving the linear system Ax = b, for x, where" );
14    let mut a = Mat64::new( 3, 3, 0.0 );
15    a[(0,0)] = 1.0; a[(0,1)] = 1.0; a[(0,2)] =   1.0;
16    a[(1,0)] = 0.0; a[(1,1)] = 2.0; a[(1,2)] =   5.0;
17    a[(2,0)] = 2.0; a[(2,1)] = 5.0; a[(2,2)] = - 1.0;
18    println!( "  * A ={}", a );
19    let b = Vec64::create( vec![ 6.0, -4.0, 27.0 ] ); 
20    println!( "  * b^T ={}", b );
21    let x = a.clone().solve_basic( &b );
22    println!( "  * gives the solution vector");
23    println!( "  * x^T ={}", x );
24    println!( "---------------------------------------------------" );
25    println!( "  * Lets solve the complex linear system Cy = d where" );
26    let mut c = Matrix::<Cmplx>::new( 2, 2, Cmplx::new( 1.0, 1.0 ) );
27    c[(0,1)] = Cmplx::new( -1.0, 0.0 );
28    c[(1,0)] = Cmplx::new( 1.0, -1.0 );
29    println!( "  * C ={}", c );
30    let mut d = Vector::<Cmplx>::new( 2, Cmplx::new( 0.0, 1.0 ) );
31    d[1] = Cmplx::new( 1.0, 0.0 );
32    println!( "  * d^T ={}", d );
33    println!( "  * Using LU decomposition we find that ");
34    let y = c.clone().solve_lu( &d );
35    println!( "  * y^T ={}", y );
36    println!( "---------------------------------------------------" );
37    println!( "  * We may also find the determinant of a matrix" );
38    println!( "  * |A| = {}", a.determinant() );
39    println!( "  * |C| = {}", c.determinant() ); 
40    println!( "  * or the inverse of a matrix" );
41    let inverse = a.inverse();
42    println!( "  * A^-1 =\n{}", inverse );
43    println!( "  * C^-1 =\n{}", c.inverse() );
44    println!( "  * We can check this by multiplication" );
45    println!( "  * I = A * A^-1 =\n{}", a * inverse );
46    println!( "  * which is sufficiently accurate for our purposes.");
47    println!( "-------------------- FINISHED ---------------------" );
48}
49
linear_systems/linear_systems.rs

linear_systems/
linear_systems.rs
