% ### LINode45EX.m ### 01.26.16 CB % Numerically integrate a general 2nd order linear autonomous system (w/ % const. coefficients) % x' = a*x + b*y % y' = c*x + d*y clear % ----------------------------------------------------- % User input (Note: All paramters are stored in a structure) P.y0(1) = 1.0; % initial value for x P.y0(2) = 1; % initial value for y P.A= [-3.9 3; % matrix A to contain coefficients A= [a b -2 1]; % c d] % Integration limits P.t0 = 0.0; % Start value P.tf = 10.0; % Finish value P.dt = 0.01; % time step % ---------------------------------------------------------------------- % +++ % determine some basic derived quantities p= P.A(1,1)+ P.A(2,2); % Tr(A) q= P.A(1,1)* P.A(2,2)-P.A(1,2)* P.A(2,1); % det(A) disp(['Tr(A)= ' num2str(p),' and det(A)= ',num2str(q)]); eigV1= [0.5*(p+sqrt(p^2-4*q)) 0.5*(p-sqrt(p^2-4*q))]; % calc. eigenvalues directly eigV2= eig(P.A); % calculate via Matlab's built-in routine disp(['eigenvalues= ' num2str(eigV1(1)),' and ',num2str(eigV1(2))]); % +++ % use built-in ode45 to solve [t y] = ode45('LINfunction', [P.t0:P.dt:P.tf],P.y0,[],P); % ------------------------------------------------------ % visualize % NOTE (re variable naming): x=y(1) and y=y(2) figure(1); clf; plot(t,y(:,1)); hold on; grid on; xlabel('t'); ylabel('x(t)') % Phase plane figure(2); clf; plot(y(:,1), y(:,2)); hold on; grid on; xlabel('x(t)'); ylabel('y(t)') % "solution space" figure(3); clf; plot(p,q,'rx','MarkerSize',9,'LineWidth',3); hold on; grid on; if (abs(p)<1), pSpan= linspace(-1,1,100); else pSpan= linspace(-1.5*p-0.1,1.5*p+0.1,100); end qSpan= pSpan.^2/4; plot(pSpan,qSpan,'k-','LineWidth',2); %ylim([-max(qSpan) max(qSpan)]) plot(pSpan,zeros(numel(pSpan),1),'b--','LineWidth',2); xlabel('Tr(A)'); ylabel('det(A)')