Trivial Example --------------- ( x^3 - y^3 ) / ( x - y ) // InputForm Cancel [ ( x^3 - y^3 ) / ( x - y ) ] // InputForm Solve [ x^3 - a + 1 == 0 , x ] // InputForm Solve [ x^5 + a * x^3 + 1 == 0 , x ] // InputForm Simple Matrix Algebra --------------------- a = { { 4.2 , 2.2 + p , -3.9 , 9.3 , 0.1 } , { 8.6 - p , r , 0.7 , -2.3 , -0.3 } , { 8.4 , -5.9 , -8.1 + q , 9.6 , 3.8 } , { -0.8 , -9.4 , -9.9 , 9.9 , 5.0 - r } , { -1.3 , -8.1 , 0.6 , -9.2 + r , -7.3 } b = { -6.8 , 2.3 , 2.7 , -7.0 , 2.0 c = { p , 1.0 , p + q , 1.0 , q a . b How Mathematica Works --------------------- FullForm [ 1.23 * x + y + 4.56 * z ] Assignment ---------- x = 1.23 * a + 4.56 y = 7.89 * x + 0.12 + z x = 1.23 * a + 4.56 x = 1.23 * x + 4.56 x =. (* Ensure that x has no value *) x = 1.23 * x + 4.56 Rules ----- rule = x ^ k_ /; k > 3 -> 0 Expand [ ( 1 + x ) ^ 7 ] /. rule rule = x ^ k_ /; k > 3 -> 0 Expand [ ( 1 + x ) ^ 7 ] /. rule Remember The Tree? ------------------ x * y + x * z /. ( y + z ) -> 0 x y + x z Delayed Rules ------------- x * y + x * z /. ( y + z ) :> 0 Worked Example -------------- x = c . MatrixPower [ a , 3 ] . c ; p = p1 * e ; q = q1 * e ; r = r1 * e ; rule = e ^ k_ /; k > 1 -> 0 y = x /. rule p = p1 * e ; q = q1 * e ; r = r1 * e ; rule = e ^ k_ /; k > 1 -> 0 y = Expand [ x ] /. rule p =. ; q =. ; r =. ; z = y /. e p1 -> p /. e q1 -> q /. e r1 -> r p = p1 * e ; q = q1 * e ; r = r1 * e ; rule = e ^ k_ /; k > 2 -> 0 y = Expand [ x ] /. rule p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y p1 =. ; q1 =. ; r1 =. ; Aside: Syntactic Gotcha ----------------------- arg = { 1 , 2 , 3 arg = { 1 , 2 , 3 fred [ ] := ( x = arg [[ 1 ]] y = x ) Using MatrixPower ----------------- For [ k = 0 , k < 21, ++k , Print [ Timing [ x = c . MatrixPower [ a , k ] . c ; p = p1 * e ; q = q1 * e ; r = r1 * e ; rule = e ^ k_ /; k > 1 -> 0 ; y = Expand [ x ] /. rule ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y ; p1 =. ; q1 =. ; r1 =. ; z ] ] ] Back To Basics -------------- PURE FUNCTION exponent ( arg ) . . . exponent = 1.0 ; value = 1.0 DO k = 1 , HUGE ( k ) value = value * arg / k temp = exponent + value IF ( temp == exponent ) EXIT END DO END FUNCTION exponent matexp [ z_ , r_ ] := Module [ { k, w, x, y } , x = z ; y = IdentityMatrix [ Dimensions [ z ] [[ 1 ]] ] ; For [ k = 1 , k < 50 , ++k , x = Expand [ x . z / k ] /. r ; w = y + x ; (* No rule needed *) y = If [ SameQ [ y, w ] , Break [ ] , w ] ] ; y ] Timing That ----------- Timing [ p = p1 * e ; q = q1 * e ; r = r1 * e ; a1 = a ; rule = e ^ k_ /; k > 1 -> 0 ; x = c . matexp [ a1, rule ] . c ; y = Expand [ x ] /. rule ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y ; p1 =. ; q1 =. ; r1 =. ; z ] Second order ------------ Timing [ p = p1 * e ; q = q1 * e ; r = r1 * e ; a1 = a ; rule = e ^ k_ /; k > 2 -> 0 ; (* Changed *) x = c . matexp [ a1, rule ] . c ; y = Expand [ x ] /. rule ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; z = y ; p1 =. ; q1 =. ; r1 =. ; z ] MatrixPower ----------- matpower [ z_ , n_ , r_ ] := Module [ { k, x, y } , k = n ; x = z ; y = IdentityMatrix [ Dimensions [ z ] [[ 1 ]] ] ; While [ k > 0 , y = If [ OddQ [ k ] , Expand [ y . x ] /. r , y ] ; k = IntegerPart [ k / 2 ] ; x = Expand [ x . x ] /. r ] ; y ] Fourier Transforms ------------------ a = { 0.906 , 4.528 , 4.693 , 1.515 + p , -2.190 , -3.239 , -0.724 , 3.238 , 5.253 - 2 * p , 3.463 , -0.906 , -4.528 , -4.693 , -1.515 + p , 2.190 , 3.239 , 0.724 , -3.238 , -5.253 - 2 * p , -3.463 Fourier [ a ] n = 20 w = Cos [ 2.0 * Pi / n ] + I * Sin [ 2.0 * Pi / n ] ; t = Table [ w ^ ( i * j ) , { i , 0 , n-1 } , { j , 0 , n-1 } ] / Sqrt [ 1.0 * n ] ; Collect [ t . a , p ] Solving Equations ----------------- a = { { 4.2 , 2.2 + p , -3.9 , 9.3 , 0.1 } , { 8.6 - p , r , 0.7 , -2.3 , -0.3 } , { 8.4 , -5.9 , -8.1 + q , 9.6 , 3.8 } , { -0.8 , -9.4 , -9.9 , 9.9 , 5.0 - r } , { -1.3 , -8.1 , 0.6 , -9.2 + r , -7.3 } c = { p , 1.0 , p + q , 1.0 , q v = LinearSolve [ a , c ] (* That's not pretty *) p = p1 * e ; q = q1 * e ; r = r1 * e ; rule1 = e ^ k_ /; k > 1 -> 0 w = ExpandAll [ v ] /. rule1 w = ExpandAll [ Together [ v ] ] /. rule1 rule2 = k1_ / ( k2_ + k3_ * e ) -> k1 * ( 1 - ( k3 / k2 ) * e ) / k2 x = w /. rule2 x = w //. rule2 y = Expand [ x ] /. rule1 p =. ; q =. ; r =. ; z = y /. e p1 -> p /. e q1 -> q /. e r1 -> r v = LinearSolve [ a , c ] ; p = p1 * e ; q = q1 * e ; r = r1 * e ; rule1 = e ^ k_ /; k > 1 -> 0 ; w = ExpandAll [ Together [ v ] ] /. rule1 ; rule2 = k1_ / ( k2_ + k3_ * e ) -> k1 * ( 1 - ( k3 / k2 ) * e ) / k2 ; x = w //. rule2 ; y = Expand [ x ] /. rule1 ; p =. ; q =. ; r =. ; z = y /. e p1 -> p /. e q1 -> q /. e r1 -> r p = p1 * e ; q = q1 * e ; r = r1 * e ; s = Expand [ a . z ] /. rule1 ; p =. ; q =. ; r =. ; t = s /. e p1 -> p /. e q1 -> q /. e r1 -> r t /. x_ /; Abs[x] < 1.0*^-12 -> 0 Second Order ------------ v = LinearSolve [ a , c ] ; p = p1 * e ; q = q1 * e ; r = r1 * e ; rule1 = e ^ k_ /; k > 2 -> 0 ; w = ExpandAll [ Together [ v ] ] /. rule1 ; rule2 = k1_ / ( k2_ + k3_ * e ) -> k1 * ( 1 - ( k3 / k2 ) * e ) / k2 ; rule3 = k1_ / ( k2_ + k3_ * ( e | e ^ * k4_ ) -> k1 * ( 1 - ( k3 / k2 ) * e - ( ( k3 / k2 ) ^ 2 ) * e ^ 2 - ( k4 / k2 ) * e ^ 2 ) / k2 ; x = Collect [ w , e ] //. rule2 //. rule3 invert [ arg_ ] := Module [ { x, y , z } , x = CoefficientList [ arg , e ] ; y = x [[ 2 ]] / x [[ 1 ]] ; z = x [[ 3 ]] / x [[ 1 ]] ; (1.0 / x [[ 1 ]] ) * ( 1.0 - y * e + ( y ^ 2 - z ) * e ^ 2 ) ] rule2 = k1_ / k2_ :> k1 * invert [ k2 ] ; x = Expand [ w //. rule2 ] /. rule1 p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; y = x ; p1 =. ; q1 =. ; r1 =. ; y p = p1 * e ; q = q1 * e ; r = r1 * e ; s = Expand [ a . y ] /. rule1 ; p =. ; q =. ; r =. ; p1 = p / e ; q1 = q / e ; r1 = r / e ; t = s ; p1 =. ; q1 =. ; r1 =. ; Print [ t ] ; t /. x_ /; Abs[x] < 1.0*^-12 -> 0 Performance ----------- For [ k = 1 , k < 10, ++ k , a = Table [ RandomReal [ ] , { k } , { k } ] + Table [ RandomReal [ ] , { k } , { k } ] * p ; b = Table [ RandomReal [ ] , { k } , { k } ] ; Print [ k , Timing [ LinearSolve [ a , b ] ; ] ] ] Eigenvalues etc. ---------------- CharacteristicPolynomial [ HankelMatrix [ 5 ] , x ] Solve [ CharacteristicPolynomial [ HankelMatrix [ 4 ] , x ] == 0 , x ] Solve [ CharacteristicPolynomial [ HankelMatrix [ 5 ] , x ] == 0 , x ] Where I Came In --------------- factorial [ n_ ] = Sqrt [ 2 * pi * n ] * ( n / e ) ^ n * (1 + 1 / ( 12 * n ) + 1 / ( 288 * n ^ 2 ) ) ; choose [ n_ , k_ ] := factorial [ k ] * factorial [ n - k ] / factorial [ n ] ; choose [ K , k ] * choose [ N - K , n - k ] / choose [ N , K ]