(* 
	AME 521-- JM Powers, Fall 1999
	
	Here we consider some stability properties of linear ODE's and
	their numerical approximation .   If we have the system
	
	dy/dt = A . y,
	
	we know that the system itself is stable if the real parts of all of the eigenvalues of A are negative .   So
	first define an arbitrary matrix A, and find its eigenvalues, λ :  
*)

A  =  {{-2, 1, -1}, {3, -4, 4}, {-1, 4, -5}};  (* Arbitrary matrix *)
A // MatrixForm
Eigenvalues[A] // N  // MatrixForm

  
\(y = {y1[t], y2[t], y3[t]}; 
f\  = \ A . y; 

    eqns\  = \ {\(y1'\)[t] == f[\([1]\)], 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        \(y2'\)[t] == f[\([2]\)], 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        \(y3'\)[t] == f[\([3]\)]}\)\(sol = DSolve[eqns, {y1, y2, y3}, t]; 
y1[t] /. sol\  // N\)y2[t] /. sol // Ny3[t] /. sol // N

h  =   .1;

(* 
	Now construct B, the matrix which is relevant to the explicit finite 
                difference technique, where we have
			
			(Subscript[y, n + 1]) = (((I + h A))^n) (Subscript[y, o]) = (B^n) (Subscript[y, o])
					
 For stability of the finite difference technique, we need for the eigenvalues of B to have modulus 
                less than one .   Note if the original problem is 
                unstable, we should expect the numerical method to be unstable 
                for sufficiently small h
	 *)

B  =  IdentityMatrix[3] + h A;
B // MatrixFormEigenvalues[B]  // MatrixForm     (* these need to have modulus  <  1 for stability *)

(*  
	Now, let' s check if the finite difference method actually 
              converges .   With n as the number of iterations, and c[[1]] as an arbitrary seed value, let' s iterate by letting B operate on c[[1]] n 
              times .   If the numerical method is stable, which requires that the real part of the eigenvalues of B have modulus less than one, the result should approach 0.  
*)
			
n  =  1000;
c  =  Table[0, {i, n}];
c[[1]]  =  {1, 5, 3};       (* seed with arbitrary value *)
Do[c[[i + 1]] = B . c[[i]], {i, n - 1}]c[[1]] // MatrixFormc[[n/2]] // MatrixFormc[[n - 1]] // MatrixFormc[[n]] // MatrixForm