. For computer algebra systems or numerical workbenches, matrix inversion is the easiest. EXAMPLE 8. Solve the initial value problem u' = Matrix([[3,-2,0],[-1,3,-2],[0,-1,3]]) u, u(0) = <0,2,6> ANSWER: u(t)=Q(t)<2,-3,1> where Q(t) =|exp(3*t)*<2,0,-1>|exp(5*t)*<2,-2,1>> METHOD: The columns of Q are solutions of the differential equation u' = Matrix([[3,-2,0],[-1,3,-2],[0,-1,3]]) u and det(Q)=-16 exp(9t) is nonzero, so Q has independent columns. The u(t)=Q(t)c for some constant vector c. We find c by solving Q(0)c=u(0)=<0,2,6>. Here's how with maple: Q0:= <<2,2,1>|<2,0,-1>|<2,-2,1>>; c:=Q0^(-1) . <0,2,6>; # Dot is matrix multiply. # c:=<2,-3,1>; EXAMPLE. Test the columns of Q(t) to see if they are solution vectors of the differential equation: Q(t) = |exp(3*t)*<2,0,-1>|exp(5*t)*<2,-2,1>> u' = Matrix([[3,-2,0],[-1,3,-2],[0,-1,3]]) u ANSWER: We test Q'(t)=P Q(t), which tests all three columns at once. It is easiest in a computer algebra system. Here's how for maple: Q:=unapply( |exp(3*t)*<2,0,-1>|exp(5*t)*<2,-2,1>>,t); LHS:=map(diff,Q(t),t); # Find the derivative of Q(t) P:=Matrix([[3,-2,0],[-1,3,-2],[0,-1,3]]); RHS:=P.Q(t); # Evaluate matrix product P times Q simplify(LHS-RHS); # The sides match if this is zero NONHOMOGENEOUS SYSTEMS We solve the vector-matrix system u'=P(t)u+f(t) where P is an nxn matrix of continuous functions on interval J. THEOREM 4. Superposition for nonhomogeneous systems The general solution of u'=P(t)u+f(t) is of the form u=uh+up, where: 1. The general solution uh(t) of u'=P(t)u is a linear combination of n independent solutions. We can write uh(t)=Q(t)c, where the columns of Q(t) are independent solutions. 2. A particular solution up(t) of u'=P(t)u+f(t) is assumed. It can be any solution whatsoever found by any method. In expanded form the general solution looks like u(t) = c1 v1(t) + ... + cn vn(t) + up(t) where c1, ..., cn are constants and v1, ..., vn are n independent solutions of the homogeneous equation u'=P(t)u. EXAMPLE 9. Consider the non-homogeneous system u'=P(t)+f(t) given by u' = Matrix([[3,-2,],[-1,3,-2],[0,-1,3]]) u + <13,-15,7> + <-9t,7t,-6t> A particular solution is xp(t) = <3t,5,2t> and the general solution of the homogeneous system u'=P(t)u is given by uh(t) = c1 exp(t) <2,2,1> + c2 exp(3t) <2,0,-1> + c3 exp(5t) <2,-2,1> The general solution in scalar form for u'=P(t)u+f(t) is x1(t) = 2 c1 exp(t) + 2 c2 exp(3t) + 2 c3 exp(5t) + 3t, x2(t) = 2 c1 exp(t) - 2 c3 exp(5t) + 5, x3(t) = c1 exp(t) - c2 exp(3t) + c3 exp(5t) + 2t, 5.2 The Eigenvalue Method for Linear Systems ==== Assumed here is basic eigenanalysis from Linear Algebra 2270. The textbook has a brief review. We study only the case of a homogeneous constant-coefficient system u'=Au where A is constant nxn. The methods apply only to diagonalizable matrices, that is, nxn matrices that have exactly n eigenpairs. A method for non-diagonalizable matrices appears in 5.5. In addition, we study work-arounds for cases not covered by section 5.3, in order to address how to solve u'=Au for all possible matrices A. THEOREM 1. The Eigenvalue Method for u'=Au Let (lambda,v) be an eigenpair of A. Then the Euler solution u(t)=exp(lambda t) v is a nonzero solution of the system u'=Au. Conversely, if the Euler equation u(t) = exp(lambda t) v provides a nonzero solution of the system u'=Au, then (lambda,v) is an eigenpair of matrix A. ALGORITHM: The Eigenvalue Method To solve u'=Au, which means "find the general solution" 1. Solve the characteristic equation |A - lambda I| = 0 for all roots lambda. 2. For each lambda in step 1, find all corresponding eigenvectors v. 3. Check that there are exactly n eigenpairs (lambda,v). If not, then this method fails to apply and no answer has been found. Otherwise, the general solution of the system u'=Au is given by u(t) = c[1] exp(lambda[1] t) v[1] + ... + c[n] exp(lambda[n] t) v[n] where (lambda[j],v[j]) are the eigenpairs and c[1],...,c[n] are constants. THEOREM. The solutions exp(lambda[j] t) v[j] are independent. REMARK. This is true even if the matrix A fails to be diagonalizable, although in this case no general solution formula was found. THEOREM. If the eigenvalues are real and distinct, then there are n eigenpairs and the general solution given in the algorithm can be assumed to have real constants in the general solution. EXAMPLE. Find the general solution by the Eigenvalue Method: x1' = 4 x1 + 2 x2, x2' = 3 x1 - x2 An equivalent formulation is u' = Matrix([[4,2],[3,-2]]) u ANSWER: u(t)=c1 exp(-2t) <1,-3> + c2 exp(5t) <2,1> WARNING: The shortcuts used in the textbook's eigenanalysis are historical and assume a background in linear algebra that has matured for a year or so. If you have no such background, then giggle at the book's shortcut and find the eigenpairs slowly by your own linear algebra skills, independent of the textbook. EXAMPLE. COMPARTMENT ANALYSIS and BRINE TANKS Let tank 1 flow into tank 2 and then into tank 3, all flow rates being r=10 gal/min. The tank volumes are V1=20, V2=40, V3=50 gallons. The initial data for the salt amounts x1, x2, x3 in pounds in the three tanks is x1(0)=15, x2(0)=0, x3(0)=0 Find the amount of salt in each brine tank for all t>0. ANSWER: u'=Au is the model where A:=Matrix([[-0.5,0,0],[0.5,-0.25,0],[0,0.25,-0.2]]) This is done by compartment analysis using dx/dt=flow in - flow out. The initial condition x1(0)=15, x2(0)=0, x3(0)=0 translates to the vector condition u(0)=<15,0,0>. ALGORITHM. We get a general solution u(t) = c1 exp(-t/2) <3,-6,5> + c2 exp(-t/4) <0,1,-5> + c3 exp(-t/5) <0,0,1> INITIAL DATA. We solve for c1, c2, c3 using u(0)=<15,0,0>, which gives the matrix equation <<3,-6,5>|<0,1,-5>|<0,0,1>> = <15,0,0> Then c1=5, c2=30, c3=125 by matrix inversion or the RREF of the augmented matrix <<3,-6,5>|<0,1,-5>|<0,0,1>|<15,0,0>> COMPLEX EIGENVALUES Technically, every formula developed earlier works for COMPLEX EIGENVALUES. The only trouble is that the solutions are complex and involve complex exponentials satisfying Euler's formula exp(at + ibt)=exp(at)(cos(bt) + i sin(bt)), i =sqrt(-1) The eigenvectors v also have COMPLEX COMPONENTS, leading to eigenvectors of the form v = c + i d, c, d = fixed vectors with real components THEOREM. To an eigenpair (a+ib,v) and its conjugate eigenpair (a-ib,w) correspond two real-valued independent solutions, constructed only from (a+ib,v) [we don't compute (a-ib,w)]: v = vReal + i vImag, where vReal, vImag are fixed vectors with real components. Solution 1 = exp(at) cos(bt) vReal - exp(at) sin(bt) vImag Solution 2 = exp(at) cos(bt) vImag + exp(at) sin(bt) vReal REMARK. These formulas replace the tedious process of taking real and imaginary parts of a complex-valued eigen-solution exp(at+ibt)(vReal+ivImag) EXAMPLE. Find the two independent solutions of u'=<<4,3>|<-3,4>>u for eigenvalue lambda=4+3i. ANSWER: The eigenpair we work with is (4+3i,<1,-i>). Then v = <1,-i> = <1,0>+i<0,-1> implies vReal = <1,0>, vImag = <0,-1> Then the solutions are sol 1 = exp(4t)cos(3t)<1,0>-exp(4t)sin(3t)<0,-1>, sol 2 = exp(4t)cos(3t)<0,-1>+exp(4t)sin(3t)<1,0> REMARK. The two solutions are formally obtained from the matrix product exp(4t) Matrix([[cos(bt),-sin(bt)],[sin(bt),cos(bt)]]) which is expanded as though vReal, vImag were scalars. The matrix is the standard rotation matrix for angle theta=bt and it represents the complex number cos(theta)+i sin(theta) as a 2x2 matrix. EXAMPLE. Re-circulating brine tanks Three brine tanks of volumes V1=50, V2=25, V3=50 gallons recirculate at rate 10 gal/min. The salt amount in pounds x1,x2,x3 in the tanks satisfies the differential system u' = Matrix([[-0.2,0,0.2],[0.2,-0.4,0],[0,0.4,-0.2]]) u where u= Find u(t) for any set of initial data. ANSWER: Let c1, c2, c3 be any constants. Then u(t) = c1 v1 + c2 v2(t) + c3 v3(t) v1 = <2,1,2> for eigenvalue lambda=0 For eigenvalue lambda=-2/5+i/5 and its conjugate v2(t) = exp(-2t/5)(cos(t/5) vReal - sin(t/5) vImag) v3(t) = exp(-2t/5)(cos(t/5) vImag + sin(t/5) vReal) vReal + i vImag = <1,-i,-1+i> = <1,0,-1>+i<0,-1,1> which implies vReal=<1,0,-1>, vImag=<0,-1,1> PHYSICAL INTERPRETATION The limits of x1, x2, x3 are the steady-state salt contents in the three tanks. MAPLE ANSWER CHECK This answer check finds a 3x3 matrix of independent columns which are solutions of the differential equation. When eigenanalysis succeeds, it returns the same answers you would get by hand, except possibly the columns of Q have been multiplied by a constant - a minor issue, they are still independent solutions. # Maple, using package DEtools and function matrixDE A:=Matrix([[-1/5,0,1/5],[1/5,-2/5,0],[0,2/5,-1/5]]); sol:=DEtools[matrixDE](A,t); Q:=Matrix(sol[1]); 5.4 Second-order mechanical systems ==== THREE FRICTIONLESS COUPLED SPRING-MASSES m1 x1'' = -k1 x1 + k2 (x2 - x1) m2 x2'' = -k2 (x2-x1) + k3 (x3 - x2) m3 x3'' = -k3 (x3 - x2) - k4 x3 M=Matrix([[m1,0,0],[0,m2,0],[0,0,m3]]) K=Matrix([[-k1-k2,k2,0],[k2,-k2-k3,k3],[0,k3,-k3-k4]]) u(t)= M u''(t) = K u(t) The vector-matrix model for 3 masses u''(t) = A u(t) Converted form, A = M^(-1) K SOLUTION of SECOND ORDER SYSTEMS THEOREM. Let (lambda,v) denote an eigenpair of nxn matrix A. Let z denote a complex number, z=a+ib. Then u(t) = exp(zt) v is a solution of u''(t) = A u(t) if and only if lambda = z^2. THEOREM. Assume matrix A is nxn and it has a real eigenvalue lambda = - omega^2, with omega positive. Then to each corresponding eigenpair (lambda,v) of matrix A corresponds two real-valued independent solutions of u'' = Au, constructed as follows: solution 1 = cos(omega t) v, solution 2 = sin(omega t) v. If lambda=0 is an eigenvalue, then the solutions are replaced by solution 1 = v, solution 2 = t v. THEOREM 1. Let nxn matrix A have distinct real negative eigenvalues lambda[1]=-omega[1]^2, ..., lambda[n]=-omega[n]^2 with corresponding eigenvectors v[1], ..., v[n]. Then the general solution of u''(t) = A u(t) is given by the expression u(t) = (a[1] cos(omega[1] t) + b[1] sin(omega[1] t)) v[1] + ... + (a[n] cos(omega[n] t) + b[n] sin(omega[n] t)) v[n] where a[1], b[1], ..., a[n], b[n] are arbitrary constants. If zero is one of the eigenvalues lambda[j], then the corresponding term in the solution is (a[j] + b[j] t) v[j] EXAMPLE. Two masses with two springs Solve Mu''=Ku for M:=Matrix([[2,0],[0,1]) and K:=Matrix([[-150,50],[50,-50]]), which corresponds to m1=2, m2=1, k1=100, k2=50, k3=0 in the standard model. ANSWER: Convert to u''=Au, A:=Matrix([[-75,25],[50,-50]]). Then apply Theorem 1 to get u(t) = (a[1] cos( 5t) + b[1] sin( 5t)) v[1] + (a[2] cos(10t) + b[2] sin(10t)) v[2] where v[1]=<1,2>, v[2]=<1,-1>. MAPLE ANSWER CHECK A:=Matrix([[-75,25],[50,-50]]); LinearAlgebra[Eigenvectors](A); # prints eigenpairs of A, # (-5^2,<1/2,1>), (-10^2,<-1,1>) EXAMPLE. Three railway cars with buffer springs that react under compression and disengage under stretch. Two cars on the right are engaged and at rest. The leftmost car travels right at velocity v0>0, hitting at time t=0 the two cars. Find the equations of motion. MODEL: u'' = Au, A:=Matrix([[-4,4,0],[6,-12,6],[0,4,-4]]) EIGENPAIRS of A: (-0^2,<1,1,1>), (-2^2,<1,0,-1>), (-4^2,<1,-3,1>) These can be found by hand analysis and checked by machine. INITIAL DATA: u(0)=<0,0,0>, u'(0)= where v0=velocity of the leftmost car (position x1, velocity x1'). DETAILS: The model holds while the buffer springs are compressed, which is while x2(t) - x1(t) < 0 and x3(t) - x2(t) < 0. Computing, they hold until t=Pi/2 and then x1=x2=x3=3 Pi v0/16, x1'=x2'=0 and x3'=v0. For t>Pi/2, speed v0 of car 1 is transferred to car 3, and cars 1, 2 are at rest. We see car 1,2 standing still and car 3 moving right with linear motion v0t+c. ANSWER: For 0 Pi/2, x1'=x2'=0 (at rest), x1=x2= (3 Pi v0)/16 and x3(t)=v0t+c with c = - 5 Pi v0/16. FORCED OSCILLATIONS and RESONANCE Consider u'' = Au + cos(omega t)F0 DEF. RESONANCE for this equation means unbounded solutions. We apply undetermined coefficients to a particular solution. THEOREM. A particular solution is up(t) = (cos(omega t))c where vector c is determined by the equation (A+omega^2I)c=-F0. The vector c is uniquely determined provided -omega^2 is not an eigenvalue of A. EXAMPLE. Two-mass system with external periodic force. Find the equations of motion and discuss resonance. MODEL: u'' = Au + (cos(omega t))F0 A:=Matrix([[-75,25],[50,-50]]), F0:=<0,50> It arises from masses m1=2, m2=1 and Hookes constants k1=100, k2=50, k3=0, with external force 50 cos(omega t) on mass m2. ANSWER: The solution considered is up(t)=(cos(omega t))c with amplitude |c|=sqrt(c1^2+c2^2) and frequency omega. The textbook uses a frequency-amplitude plot to detect resonance at omega^2 equal to 25 or 100, which are negatives of the two eigenvalues of matrix A. For the cases omega=5 and omega=10, the solutions are unbounded because of the term t sin(omega t) in the answer. For omega=5: x1(t)= (5/3)(t)sin(5t)-(7/18)cos(5t) x2(t)=(10/3)(t)sin(5t)- (1/9)cos(5t) For omega=10: x1(t)=-(67/72)cos(10t)-(5/6)(t)sin(10t) x2(t)= (19/72)cos(10t)+(5/6)(t)sin(10t) MAPLE STUDY We employ dsolve() to find the solution when omega^2=25 or 100. The result in each case is an unbounded solution. A:=Matrix([[-75,25],[50,-50]]); B:=cos(omega*t)*<0,50>; u:= ; sys:=A.u+B-map(diff,u,t,t); dsolve([sys[1]=0,sys[2]=0],[x1(t),x2(t)]); sys1:=subs(omega=5,sys); dsolve([sys1[1]=0,sys1[2]=0],[x1(t),x2(t)]); sys2:=subs(omega=10,sys); dsolve([sys2[1]=0,sys2[2]=0],[x1(t),x2(t)]); PERIODIC and TRANSIENT SOLUTIONS DEF. The TRANSIENT SOLUTION is the set of terms in the general solution which damp to zero at t=infinity. They are terms in the homogeneous solution for mechanical systems u''=Au+F0 cos(omega t) DEF. The STEADY-STATE SOLUTION is the set of terms remaining in the general solution of u''=Au+F0 cos(omega t) after all terms have been removed that limit to zero at t=infinity. DEF. The STEADY PERIODIC SOLUTION is the steady-state solution in the special case when it is T-periodic. THEOREM. For the mechanical system u'' = Au + F0 cos(omega t) where matrix A has only negative eigenvalues, there is a unique steady periodic solution xp(t) obtained by undetermined coefficients, for -omega^2 not an eigenvalue of A. EARTHQUAKES MODEL: u'' = Au + (E omega^2 cos(omega t)) b A=M^(-1) K M=diagonal matrix of masses [see equation (5) in 5.4] K=Hooke's matrix [see equation (6) in 5.4] b=<1,1,...,1> E=amplitude of the ground oscillation omega=frequency of the ground oscillation MAPLE STUDY The problem can be studied for any number of floors in the building, using the 2x2 example recorded above. ===== end