% problem notes Ch8, S2009

%24 Apr, 8.1: 4, 12, 38 
%24 Apr, 8.2: 4, 19 

=======
8.1-4: 
=======
    [ 3  -1]
A = [      ]
    [ 1   1]
Cayley-Hamilton-Ziebur method: The eigenvalues are 2,2 and therefore

        x(t) = c1 exp(2t) + c2 t exp(2t)

That is, Euler's theorem with a double root r=2 implies two atoms
exp(2t), t exp(2t) and component x(t) must be a linear combination of these two
atoms. 

The first DE x'=3x-y implies 
     y = 3x-x'
       = 3c1 exp(2t)+3c2 t exp(2t) - 2c1 exp(2t)
         -c2 exp(2t) - 2c2 t exp(2t)= (c1-c2) exp(2t) + c2 t exp(2t). 
Then
                x(t) = c1 exp(2t) + c2 t exp(2t)
                y(t) = (c1-c2) exp(2t) + c2 t exp(2t)

Take the partials on symbols c1, c2 to find a fundamental matrix

              [ exp(2t)  t exp(2t)        ]
        Phi = [                           ]
              [ exp(2t)  t exp(2t)-exp(2t)]

The general solution of u'=Au is then u=(Phi)c where c is an arbitrary
2-vector. Solve u(0)=Phi(0)c for c and then report the answer.


Laplace resolvent method: Start with u'=Au, then

                                      [ s-3  1  ][1]
      L(u)=inverse(sI - a)u(0)=inverse[         ][ ]
                                      [ -1   s-1][0]
     
            [(s-1)/(s^2-4*s+4)] [L(t*exp(2*t)+exp(2*t))]
          = [                 ]=[                      ]
            [1/(s^2-4*s+4)    ] [L(t*exp(2*t))         ]

Then 
                        x(t) = t exp(2t)+exp(2t)
                        y(t) = t exp(2t)

The fundamental matrix can be taken to be Phi=invlaplace(B) where

                                  [ s-3  1  ]
        B=inverse(sI - a) =inverse[         ]
                                  [ -1   s-1]

           [(s-1)/(s^2-4*s+4)  -1/(s^2-4*s+4)   ]
         = [                                    ]
           [1/(s^2-4*s+4)      (s-3)/(s^2-4*s+4)]

Then

              [t exp(2t)+exp(2t)  -t exp(2t)         ] 
        Phi = [                                      ] = exp(At) = Exp Matrix
              [t exp(2*t)         -t exp(2*t)+exp(2t)]

Some maple help:

 with(linalg):
 with(inttrans):

 # Find eigenvalues of A, then use Cayley-Hamilton method
 A:=matrix([[3,-1],[1,1]]);
 eigenvals(A);

 # Laplace resolvent method to find u(t)=exp(At)u(0)
 B:=inverse(evalm(s*diag(1,1)-A));
 evalm(B&*vector([1,0]));map(invlaplace,%,s,t);

 # Find exp(At)
 inverse(evalm(s*diag(1,1)-A)); map(invlaplace,%,s,t);
   # alternate is
 exponential(A,t);

 
=======
8.1-12:
=======
    [ 5  -4]
A = [      ]
    [ 3  -2]

The methods of 8.1-4 apply. In particular:

=======
METHOD 1: Solve by Laplace theory.
=======
 # Find exp(At)=inverse laplace of the resolvent. Def: Resolvent == inverse(sI -A).
 with(linalg): with(inttrans):
 A:=matrix([[5,-4],[3,-2]]);
 inverse(evalm(s*diag(1,1)-A)); Phi:=map(invlaplace,%,s,t);

 Then

                [-3*exp(t)+4*exp(2*t), 4*exp(t)-4*exp(2*t) ]
      exp(At) = [                                          ]
                [-3*exp(t)+3*exp(2*t), 4*exp(t)-3*exp(2*t) ]

=======
METHOD 2. Putzer Algorithm for exp(At)
=======
                               exp(lambda_1 t) - exp(lambda_2 t)
 exp(At) = exp(lambda_1 t) I + --------------------------------- (A - lambda_1 I)
                                      lambda_1 - lambda_2

The above is for 2x2 matrices A only. There is a general Putzer algorithm [see the manuscript on systems].

DOUBLE ROOT. When lambda_1 = lambda_2, then limit on the fraction with L'Hopital's Rule, as lambda2 --> lambda1, 
to get the replacement term t exp(lambda_1 t). 

COMPLEX ROOTS. When lambda_1 and lambda_2 are complex conjugates, then take the REAL PART of the formula, 
because exp(At) is REAL [discard all terms containing complex unit i].

For this problem, lambda1 = 1, lambda2 = 2,

                    [ 1   0 ]                       [ [ 5  -4 ]      [ 1  0 ] ]
  exp(At) = exp(t)  [       ]  +   (exp(2t)-exp(t)) [ [       ]  - 1 [      ] ]
                    [ 0   1 ]                       [ [ 3  -2 ]      [ 0  1 ] ]
=======
METHOD 3. Fundamental matrix formula exp(A t) = Z(t) inverse(Z(0))
=======
Find two independent solutions of u'=Au. Define Z(t)=augmented matrix of the
two solutions. Then do a matrix multiply to find exp(A t).

Ziebur's method, applied as the book's shortcut, finds the solution of u'=Au,
   x'(t) = 5 x - 4 y,
   y'(t) = 3 x - 2 y,
by determining the char equation roots 1,2 and then atoms exp(t), exp(2t). Then
the equation is solved by the equations
   x(t) = c1 (atom 1) + c2 (atom 2) = c1 exp(t) + c2 exp(2t),
   y(t) = (1/4)(5x - x') = c1 exp(t) + (3/4)c2 exp(2t).
The last equation is found from x'=5x-4y by solving it for y.
Then Z(t) is the matrix whose first column is the partial of this 
answer on c1, and whose second column is the partial of this answer on c2.
            [ exp(t)     exp(2t)     ]
    Z(t) =  [                        ]
            [ exp(t)    (3/4)exp(2t) ]
Finally, exp(At) = Z(t) inverse(Z(0)) with 
            [ 1     1 ]
    Z(0) =  [         ]
            [ 1   3/4 ]
gives after matrix multiplication the same answer reported above.

=======
8.1-38: 
=======
The methods of 8.1-12 apply. In particular:

=======
METHOD 1. Maple resolvent method, exp(At) = inverse laplace of the resolvent.
=======
 # Find exp(At) from the resolvent == inverse(sI - A)
 with(linalg): with(inttrans):
 A:=matrix([[5,20,30],[0,10,20],[0,0,5]]);
 inverse(evalm(s*diag(1,1,1)-A)); Phi:=map(invlaplace,%,s,t);

 Then

    [exp(5*t) 4*exp(10*t)-4*exp(5*t) -50*t*exp(5*t)-16*exp(5*t)+16*exp(10*t)]
Phi=[0        exp(10*t)               4*exp(10*t)-4*exp(5*t)                ]
    [0        0                       exp(5*t)                              ]

=======
METHOD 2. Fundamental matrix method, exp(At)=Phi(t)inverse(Phi(0)). 
=======
The matrix Phi in this case is the augmented matrix of eigenvector solutions 
exp(lambda_k t) v_k which appear Fourier's differential system theorem 
[THEOREM 3 in the book, chapter 7]:

        u(t) =c1 exp(r1 t) v1 + c2 exp(r2 t) v2 + c3 exp(r3 t) v3.

Then Phi is the agumented matrix of exp(r1 t) v1, exp(r2 t) v2, exp(r3 t) v3. 
Unfortunately, matrix A is not diagonalizable and there are not 3 eigenpairs.
So the method fails. Then the book suggests to find generalized eigenvectors.
This method goes way beyond the level of the 2250 course, so we abort the idea
and try something else.

Engineers should learn the Laplace Resolvent method and
Putzer's method, because the eigenanalysis method of Theorem 3 is
limited in normal use to matrices A with n distinct real eigenvalues.

=======
METHOD 3. Fundamental matrix, exp(At) = Z(t) inverse(Z(0))
=======
The problem reduces to finding the 3 independent columns of Z(t), which are
independent solutions of u'=Au. The equation u'=Au can be written in scalar
form as the triangular system
  x' = 5x + 20y + 0x,
  y' =      10y + 20z,
  z' =             5z.
This system can be solved by the linear integrating factor method of chapter 1.
First, verify that 
    z = c3 exp(5t) 
      = constant divided by the integrating factor.
The integratin factor is the exponential of the integral of -5. Stuff the answer 
for z into the second DE, then solve to get y(t). It will contain a constant c2, 
as well as c3. Stuff z(t), y(t) into the first DE and solve. Then x(t) contains 
a constant c1, plus constants c2, c3.

The columns of Z(t) are the partials of the scalar answer x(t), y(t), z(t) on
the symbols c1, c2, c3, respectively. This is because
     [x(t)]        [c1]
     [y(t)] = Z(t) [c2]
     [z(t)]        [c3]
Now compute Z(0), finds its inverse, and perform a matrix multipy to get the
exponential matrix answer
     exp(At)=Z(t) inverse(Z(0)).

=======
ANSWER CHECK: There is a direct route to finding exp(At) in maple:
=======
 with(linalg):A:=matrix([[5,20,30],[0,10,20],[0,0,5]]);
 exponential(A,t);

=======
8.2-4: 
=======
    [ 4  1]       [ exp(t)]         [1]
A = [     ]   F = [       ]  u(0) = [ ]
    [ 6 -1]       [-exp(t)]         [1]

FIRST SOLUTION.
Superposition says u=uh+up. 

Find exp(At) for example by using the Cayley-Hamilton-Ziebur method for Phi
followed by exp(At)=Phi(t)inverse(Phi(0)). See also the MAPLE SOLUTION below.

               [exp(-2t)+6 exp(5t)     -exp(-2t)+exp(5t)  ]
 exp(At)= (1/7)[                                          ]
               [-6 exp(-2t)+6 exp(5t)   6 exp(-2t)+exp(5t)]

Then uh=exp(At)u0, with u0 unknown just now, to be found later.

Undetermined coefficients suggests to find up=exp(t)c where c is a
2-vector with component c1, c2. Substitution into u'=Au+F gives

                         exp(t)c = exp(t)Ac + F

Then                                           [ 1]
                               c = Ac + b, b = [  ]
                                               [-1]

has unique solution c1=-1/12, c2=-3/4 The result is 

                                    [-1/12]
                         up = exp(t)[     ]
                                    [-3/4 ]

Add u=uh+up and determine u0 from the formula

                   u0 = uh(0)+up(0)
                   u0 = u(0)+vector([-1/12,-3/4]) = vector([11/12,-7/4])

Then u=uh+up=exp(At)u0+up gives
     
               [33/28*exp(5*t)-1/12*exp(t)-2/21*exp(-2*t)]
           u = [                                         ]
               [33/28*exp(5*t)-3/4*exp(t)+4/7*exp(-2*t)  ]


MAPLE SOLUTION.
We will use maple assist to find exp(At), then proceed in maple to
do the vector integrations required in variation of parameters.

 # Find exp(At)
 with(linalg): with(inttrans):
 A:=matrix([[4,1],[6,-1]]);
 inverse(evalm(s*diag(1,1)-A)); Phi:=map(invlaplace,%,s,t);

 # Integrate var of param formula up=int(G,t=0..x);
 F:=vector([exp(t),-exp(t)]); u0:=vector([1,1]);
 G:=evalm(inverse(Phi)&*F);
 H:=subs(x=t,map(int,G,t=0..x));
 up:=evalm(Phi&*H);

 # Find uh from the initial condition
 uh:=evalm(Phi&*u0);

 # Superposition u=uh+up
 u:=evalm(uh+up);
 map(simplify,%);

     [33/28*exp(5*t)-1/12*exp(t)-2/21*exp(-2*t)]
 u = [                                         ]
     [33/28*exp(5*t)-3/4*exp(t)+4/7*exp(-2*t)  ]

=======
8.2-19:
=======
Solve u'=Au+F with initial condition.

    [1  2]       [180t]         [0]
A = [    ]   F = [    ]  u(0) = [ ]
    [2 -2]       [90  ]         [0]

=======
METHOD 1. Laplace resolvent exp(At) = inverse laplace of the resolvent matrix.
=======
 # Find exp(At) by Laplace resolvent == inverse(sI -A)
 with(linalg): with(inttrans):
 A:=matrix([[1,2],[2,-2]]);
 inverse(evalm(s*diag(1,1)-A)); Phi:=map(invlaplace,%,s,t);

 =======
METHOD 2. Variation of parameters.
=======
We will use maple assist to find exp(At), then proceed in maple to do
the vector integrations required in variation of parameters. The code is
modified from the solution of the previous problem.

# Integrate var of param formula up=int(G,t=0..x);
 F:=vector([180*t,90]); u0:=vector([0,0]);
Phi:=exponential(A,t);
 Integrand:=evalm(inverse(Phi)&*F);
 Integral:=subs(x=t,map(int,Integrand,t=0..x));
 up:=evalm(Phi&*Integral);

 # Find uh from the initial condition, uh = exp(At) u0
 # uh:=evalm(Phi&*u0);
 uh:=vector([0,0]); # because of zero ic

 # Superposition u=uh+up
 u:=evalm(uh+up);
 map(simplify,%);

     [-60*t-70+16*exp(-3*t)+54*exp(2*t)]
 u = [                                 ]
     [-60*t+5-32*exp(-3*t)+27*exp(2*t) ]

%%===end of Ch8 problem notes======