% 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.