Chapter 5 Problem Notes Updated 2015 5.1-16: ========================================================================= Define u1=x(t), u2=y(t), u3=z(t). Then u' = Au for the 3x3 matrix [2 -3 0] A = [1 1 2] [0 5 -7] 5.1-24: ========================================================================= The general solution of the system is x = c1 x1 + c2 x2. Independence of x1, x2 is determined by the value of the wronskian of x1 and x2, which is the determinant of the matrix [ exp(3t) 2 exp(-2t)] W = [ ] [ -exp(3t) -2 exp(2t) ] Evaluate det(W) and show it is nonzero. The two requested answer checks are for the two differential equations x1' = Ax1 and x2' = Ax2 where x1, x2 are the columns of matrix W and [ 4 1] A = [ ]. [-2 1] The two answer checks can be done at once by verifying W' = AW. The latter is done by expanding the LHS and the RHS, showing LHS=RHS. ========================================================================= 5.2-10: ========================================================================= The system can be written u'=Au where [x1] [-3 -2] u =[ ], A = [ ] [x2] [ 9 3] The characteristic equation of A is r^2 + 9 = 0 with roots 3i, -3i. The corresponding eigenvectors are [i-1] [-i-1] v1 =[ ], v2 = [ ] [ 3 ] [ 3 ] These should be found from one frame sequence starting with B = A - (3i)I and ending with rref(B), by applying the last frame algorithm. The second eigenvector v2 is found from the THEOREM: v2 = conjugate(v1), which applies to conjugate eigenvalues only. The solution by the eigenvalue method [Edwards Penney naming convention] is then u = c1 v1 exp(3ti) + c2 v2 exp(-3ti). This formula is expected to be translated to its real form as in Example 3, section 5.2 of the textbook. The method is to focus on one complex solution v1 exp(3ti) = v1(cos 3t + i sin 3t) and then take real and imaginary parts to get two vector solutions w1, w2 which are independent. Then a REAL solution is the formula u = c1 w1 + c2 w2. Your problem is to find formulas for vectors w1 and w2 from the identity v1 exp(3ti) = v1(cos 3t + i sin 3t) = w1 + i w2. ============ Cayley-Hamilton method ========= The system can be written u'=Au where [x1] [-3 -2] u =[ ], A = [ ] [x2] [ 9 3] The characteristic equation of A is r^2 + 9 = 0 with roots 3i, -3i. Then Euler's theorem says that the atoms are cos 3t, sin 3t. The CH Method implies there are vectors c1, c2 such that u= cos 3t c1 + sin 3t c2. The remaining details find vectors c1, c2, which are not arbitrary, but uniquely determined by matrix A and vector u(0). ADVICE: Differentiate the above relation n-1 times and set t=0 to get n vector equations for the vector unknowns c1, c2, c3, c4, ... and then solve for these vector unknowns. ADVICE for this EXAMPLE: Differentiate 1 time, set t=0 to get 2 equations for the vectors c1, c2. Then solve for c1, c2. u(0)=(1) c1 + (0) c2 Au(0)=u'(0)=3(0) c1 + 3(1) c2 Assign symbols to u(0) and then determine Au(0) by matrix multiply: [a] [-3 -2][a] [-3a-2b] u(0)=[ ] ==> Au(0) = [ ][ ] = [ ] [b] [ 9 3][b] [ 9a+3b] Then the above equations for vectors c1, c2 are solved for c1, c2 to give the general solution [a] sin 3t [-3a-2b] u(t)= (cos 3t)[ ] + ------ [ ] [b] 3 [ 9a+3b] ============ Laplace Resolvent method ================= LAPLACE THEORY. This problem can also be solved by Laplace theory applied to the matrix differential system u' = Au, u(0)=u0 Step 1. Transform the system to (sI-A) L(u) = u0 Step 2. Use the resolvent of s I - A, which is inverse(sI-A), to write L(x) = inverse(s I - A)(u0) [-3 -2] A = [ ] [ 9 3] char polynomial = delta=det(sI-A)=s^2+9 1 [s-3 -2 ] resolvent = ----- [ ] delta [ 9 s+3] Then the resolvent formula for L(x) gives 1 [s-3 -2 ][a] [a] L(x) = ----- [ ][ ] where u0 = [ ] delta [ 9 s+3][b] [b] The rest of the problem is easily done by hand, preparing the 4 matrix entries for the backward Laplace table. The answer can be compared to others above, but in general, don't expect symbols a,b to have the same meaning across the various solution methods. These symbols are arbitrary constants, so it may require some work to see the answers are equivalent. ========================================================================= 5.2-16: ========================================================================= The system can be written u'=Au where [x1] [-50 20] u =[ ], A = [ ] [x2] [100 -60] The characteristic equation of A is r^2 + 110r + 1000 = 0 with roots -10, -100. The corresponding eigenvectors are [ 1 ] [ 2 ] v1 =[ ], v2 = [ ] [ 2 ] [ -5 ] These should be found from two frame sequences starting with B = A - (lambda)I [using lambda = -10 then lambda = -100] and ending with rref(B), by applying the last frame algorithm. The solution by the eigenvalue method [Edwards Penney naming convention] is then u = c1 v1 exp(-10t) + c2 v2 exp(-100t). ========================================================================= 5.3-6: ========================================================================= This problem follows exactly Example 1 in the textbook, section 5.3. An answer is provided in the answer section of the textbook. The second order matrix differential system is [1 0] [-(2+4) 4 ] [-6 4] M x'' = K x, M = [ ], K = [ ] = [ ] [0 2] [ 4 -(4+4)] [ 4 -8] Then x'' = Ax with A=inverse(M)K: [-6 4] A = [ ] [ 2 -4] The eigenpairs of A are found from two frame sequences. Show the details. [-2] [1] (-8, [ ]), (-2, [ ]) [ 1] [1] The NATURAL FREQUENCIES are sqrt(-lambda): omega1 = sqrt(8), omega2 = sqrt(2) The NATURAL MODES x1(t) and x2(t) are the two vector terms in (14), Example 1 of section 5.3, with changes made for this problem. Describe the two modes with the same language as used in the textbook. ================ LAPLACE THEORY. This problem can also be solved by Laplace theory applied to the matrix differential system x'' = Ax, x(0)=c2, x'(0)=c1. However, the solution x(t) is not in a form which identifies the natural modes. What is useful is the Laplace model itself, which identifies the natural frequencies as the roots of the denominator delta [see below] in all the Laplace entries of the resolvent. Step 1. Transform the system to s^2 L(x) = A L(x)+ c1 + c2 s. Step 2. Use the resolvent of s^2 I - A, which is inverse(s^2I-A), to write L(x) = inverse(s^2 I - A)(c1 + c2 s) delta = s^4+10*s^2+16 = (s^2+8)(s^2+2) 1 [s^2+4 4 ] resolvent = ----- [ ] delta [ 2 s^2+6] Then the resolvent formula for L(x) gives 1 [s^2+4 4 ] 1 [s^s+4s 4s ] L(x) = c1 ----- [ ] + c2 ----- [ ] delta [ 2 s^2+6] delta [ 2s s^3+6s] The rest of the problem is partial fractions, which can be done with maple assist using the function convert(). ========================================================================= 5.3-20 ========================================================================= 7.3-30: ========================================================================= Start with k1=r/V1=10/25, k2=r/V2=10/40. Then write the system in scalar form and form the coefficient matrix [-2/5 1/4] A = [ ] A = [2/5 -1/4] The characteristic equation is r^2 + 10*r/20 = 0 with roots r=0,-13/20. The eigenvectors corresponding are v1 = [5, 8], v2 = [1, -1] Initial data x1(0)=15, x2(0)=0 is applied to the eigenanalysis solution of u' = Au, u = [x1(t), x2(t)], to obtain the vector equation u(t) = c1*v1*exp(0t) + c2*v2*exp(-13t/20) Then c1=25/13, c2=120/13. The x1-scene and x2-scene graphics are exponential-like graphs, the first decreasing and the second increasing, with corresponding aymptotic limits x1=75/13, x2=120/13. ========================================================================= 5.4-24, 4th Edition ========================================================================= The matrix A should be [ -4 4 0] A = [ 16 -32 16] [ 0 4 -4] The time reported in the problem is t=Pi/2. The answers x1'(t)=-v0/9 and x2'(t)=x3'(t)=8 v0/9 are correct. Follow example 2, section 5.4. The rebound of car 1 to the left, and energy is given to cars 2,3 to move them right. The problem is analyzed with the same ideas as in example 2. Include a re-draw of figure 5.4.5. # Maple code for Exercise 5.4-24, 4th edition Edwards-Penney with(LinearAlgebra): # units == fps, but given tons M:=(2000/32)*<32,0,0|0,8,0|0,0,32>; K:=2000*<-4,4,0|4,-8,4|0,4,-4>^+; A:=M^(-1).K; Eigenvectors(A); # Solving u'' = A u with u=, u(0)=<0,0,0>, u'(0)=. # Reference: Example 2 in Edwards-Penney section 5.4 # Get a1=a2=a3=0 and then solve for b1,b2,b3 as in Example 2, 5.4 f0:= unapply(b[1]*t*<1,1,1>+b[2]*sin(2*t)*<-1,0,1>+b[3]*sin(6*t)*<1,-8,1>,t); f0(t); eq:={diff(f0(t)[1],t)=v0,diff(f0(t)[2],t)=0,diff(f0(t)[3],t)=0}; eq1:=simplify(subs(t=0,eq)); solve(eq1,{b[1],b[2],b[3]}); # left off factor v0, for simplicity and graphing. Equivalent to v0 = 1 X:= unapply((4/9)*t*<1,1,1>+(-1/4)*sin(2*t)*<-1,0,1>+(1/108)*sin(6*t)*<1,-8,1>,t); f1:=unapply(X(t)[1],t);f2:=unapply(X(t)[2],t);f3:=unapply(X(t)[3],t); # Look for t-interval on which x2-x1 < 0 and x3 - x2 < 0. # Get 0 < t < Pi/2. plot({f2(t)-f1(t)},t=0..Pi); plot({f3(t)-f2(t)},t=0..Pi); fsolve(f2(t)=f1(t),t=1..2);fsolve(f3(t)=f2(t),t=1..2); evalf(Pi/2); # So t=Pi/2 is the transition time to linear motion # Find out velocities at t=Pi/2. xprime:=evalf(subs(t=Pi/2,diff(f1(t),t))); yprime:=evalf(subs(t=Pi/2,diff(f2(t),t))); zprime:=evalf(subs(t=Pi/2,diff(f3(t),t))); # Compare answers for position and velocity with 1/9 and 8/9 evalf(1/9); evalf(8/9); # Decimals to fractions # Find out positions at t=Pi/2 xpos:=evalf(subs(t=Pi/2,f1(t)));ypos:=evalf(subs(t=Pi/2,f2(t))); zpos:=evalf(subs(t=Pi/2,f3(t))); # Answer check for u'' = Au, u(0)=<0,0,0>, u'(0)=<1,0,0> # Define u(t), u'(t), u''(t) u:=X(t); up:=: upp:=; upp-A.u; # Should be zero subs(t=0,u); # Expect zero subs(t=0,up); ========================================================================= 5.3-24 Two-axle automobile ========================================================================= 100 x'' = -(4000) x + (2000)(2) theta 1000 theta'' = (2000)(2) x - (2000)(36+16) theta The two natural frequencies of oscillation omega1, omega2 are those frequencies which appear in the atoms constructed from the roots of the characteristic equation det(A - r^2 I)=0. The symbols are defined by A = M^{-1} K K = Hooke's matrix == matrix([[-4000,2000*2],[2000*2,-2000*(6^2+4^2)]]); M = mass matrix == matrix([[100,0],[0,1000]]); There are four values for r, because the characteristic equation is a quartic. The values are in complex conjugate pairs. Euler's theorem implies there are 4 atoms, which are base trig atoms, sine and cosine. The frequencies are 6.13, 10.31, approximately. Critical speeds v1, v2 are found by solving the equations omega1 = 2 Pi v1/40 omega2 = 2 Pi v2/40 Maple Answer Check: K:= matrix([[-4000,2000*2],[2000*2,-2000*(6^2+4^2)]]); M:=matrix([[100,0],[0,1000]]); with(linalg): A:=evalm(inverse(M)&*K); solve(det(A-r^2*diag(1,1))=0,r); ========================================================================= 5.4-11 Find the general solution of u'=Au ========================================================================= Matrix A is given by -3 0 4 -1 -1 -1 1 0 1 No notes yet. 5.4-29 Find the general solution of u'=Au ========================================================================= Matrix A is given by -1 1 1 -2 7 -4 -6 11 5 -1 1 3 6 -2 -1 6 No notes yet. ========================================================================= 5.5-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); ========================================================================= 5.5-12: ========================================================================= [ 5 -4] A = [ ] [ 3 -2] The methods of 5.5-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. ========================================================================= 5.5-38: ========================================================================= The methods of 5.5-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); ========================================================================= 5.6-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) ] ========================================================================= 5.6-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) ] =========================================================================