% Using eigenvalues and vectors to solve a system of ODEs A = [0 -1; 1 0]; % matrix that defines the linear ODEs x0 = [1 0]'; % initial condition t=3*pi/2; % set the time at which you want to evaluate x [V,D] = eig(A); % find eigenvalues and vectors of A % The solution is given by the command: x_of_t = V * diag( exp( diag(D*t) ) ) * ( V \ x0 ) % ( V \ x0 ) solves for the arbitrary solution constants in terms of ICs. % diag( exp( diag(D*t) ) ) gives exp(lambda*t) on the diagonals % and zero elsewhere. % The whole command thus gives a linear combination of the eigenvectors % of A with each eigenvector multiplied by the appropriate constant % to satisfy the ICs and by e raised to the power of (the associated % eigenvalue lambda)*t Note that the solution is real but several of % the intermediate calculations are complex. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % NEWS FLASH (12/4/96): The command mistakenly used in both Ruina % and Terrell's lectures as well as in Nagy's solutions to % homework 8 (which was mistitled 'Solution to HW#7') was: x_of_t = V * exp( D*t ) * ( V \ x0 ) % THIS IS WRONG, FIX YOUR NOTES % Most students (trusting their professors as they probably shouldn't!) % used the wrong command above in their HW and got almost the right answer. % The problem is this: exp( D*t ) takes every element of D*t and raises % it to the e power. On the diagonal this gives what is desired. But % off the diagonal it gives exp(0) which is 1 (one) in every position % where what we want is zero. % Why does the wrong command give almost the right answer? % Because in the assigned problem the diagonal elements of % exp(D*t) were so big that the ones off the diagonal were effectively % almost like the desired zeros. % A command, almost looking like the wrong command above, but that happens % to be right, is: x_of_t = V * expm( D*t ) * ( V \ x0 ) %THIS IS CORRECT % Here the command EXPM is used. It is the matrix exponential, something % you can read about in EDWARDS and PENNEY section 5.5 (which we didn't % cover). The matrix exponential is NOT the matrix each of whose % components are raised to the e power. In the case of a diagonal % matrix (like D or D*t) the matrix exponential IS the diagonal % components raised to the e power with the zero's off the diagonal % being left alone. % Why did Andy write that wrong command and then convince Professor % Terrell, TA Tamas, and hundreds of students to do something wrong? % Because he wrote out the correct command with DIAG above (thats % why he put DIAG in the commands in the MATLAB matrix examples file) % and then, on a later day under the influence of a stray cosmic % ray, 'realized' (incorrectly) that exp(0)=0. This led to the % erronious conclusion that the two DIAG commands could be eliminated.