sol_ex21_and_ss15

sol2_2017

Ex22.m Objective-C source code, 2 KB

File contents

function Ex22(x0)
	% octave code for solving simultaneously the time-varying ode
	% and the matrix differential equation for the fundamental matrix
	% as well as the monodromy matrix, Floquet multipliers and Lyapunov
	% exponents.
	% Call with Ex22(x0)
	% where x0 is the initial condition for the simulation.


	% Simulation time vector
	T = 2*pi % 1 period of the matrix A
	t = linspace(0,T,200);

	% Initial condition for the fundamental matrix in vector form,
	% chosen as identity for the purpose of calculating the monodromy
	% matrix and the Floquet exponents
	p0 = [1;0;0;1];
    
	% Solving the differential equations
	[ts, sol] = ode45(@dyn,t,[x0;p0]);

	% Determination of the monodromy matrix if the simulation time T
	% has been chosen as one period!
	M = [sol(end,3), sol(end,4);
	     sol(end,5), sol(end,6)];
	[v,lamb] = eig(M);
	mu=[lamb(1,1);lamb(2,2)];
	
	% Determining the logarithm of M (Notice, that M is diagonalizable
	% with diagonal form M1..)

	M1 = inv(v)*M*v;            % Diagonalization of M
	B=zeros(2,2);               % Initialization of B = 1/T*log(M)
	B(1,1) = 1/T*log(M1(1,1));  
	B(2,2) = 1/T*log(M1(2,2));
	sigma = real(eig(B));       % Lyapunov exponents

	fprintf('Monodromy matrix:\n');
	M
	fprintf('Floquet multipliers:\n');
	mu
	fprintf('log(M):\n');
	B
	fprintf('Lyapunov exponents:\n');
	sigma

	% Graphical output
	% Time responses of the state variables
	figure(1)
	subplot(2,1,1)
	plot(ts,sol(:,1))
	ylabel('x_1')
	grid on
	subplot(2,1,2)
	plot(ts,sol(:,2))
	ylabel('x_2')
	xlabel('Time t')
	grid on
    	% Phase plane plot
	figure(2)
	plot(sol(:,1),sol(:,2))
	grid on
	hold on

	

end

function dxdt = dyn(t,x)
	u=0;% place holder for control problem..

	% Changing back to the matrix notation for P
	P=[x(3),x(4); 
	   x(5),x(6)];
	% the system dynamics matrix
	A = [-1+cos(t), 0; 1, -.5+sin(t)];
	% Determining the right-hand side of the ode 
	dxdt(1:2,1) = A*x(1:2)+[1;0]*u;	

	% and the matrix differential equation:
	dPdt = A*P;

	% completing the right-hand-side in vector notation
	dxdt(3,1) = dPdt(1,1);
	dxdt(4,1) = dPdt(1,2);
	dxdt(5,1) = dPdt(2,1);
	dxdt(6,1) = dPdt(2,2);

end

solEx54

Ex5_4.m Objective-C source code, 2 KB

File contents

function FloquetVanDerPol
clc
fprintf('Determination of the Floquet and Lyapunov exponents for the stable\n')
fprintf('limit cycle of the van der Pol oscillator\n')
fprintf('       ddt x + mu*(x^2-1)*dt x + x = 0\n')
fprintf('with mu = 2.\n')
fprintf('\n')
% ------------------
% A. Schaum, SS-2016
% ------------------
close all
% Solution of the van der Pol oscillator dynamics
t = linspace(0,50,200);
x0 = [0.1,0.2];
p.mu = 2;
[tv,solv] = ode45(@(t,x)vdp(t,x,p),t,x0);
figure(1)
plot(tv,solv(:,1),'LineWidth',2);

% Graphically determined starting and end time of one cylce
fprintf('Starting point of evaluation:\n')
t0 = 30.90455 
t1 = 38.4425;
%and period 
fprintf('Approximate period:\n')
T  = t1-t0
% Interpolation of numerical solution over the interval [t0, t0+T]
p.tf = linspace(t0,t0+T,200);
% Graphical comparison of the interpolation with the original solution
p.gam = interp1(tv,solv,p.tf);
figure(1);hold on;
plot(p.tf,p.gam(:,1),'r','LineWidth',2)

% Solution of the fundamental equation dt P=AP with A=df/dx evaluated at
% the approximated limit cycle p.gam.
% The matrix P = [p1,p2;p3,p4] is handled using the vector of its entries
% xp = [p1,p2,p3,p4]'. As starting point the identity matrix is chosen.
P0 = eye(2,2);
xp0 = [P0(1,1),P0(1,2),P0(2,1),P0(2,2)]';

[tp,solp] = ode45(@(t,x)fundeq(t,x,p),p.tf,xp0);

% The monodromy matrix is given by
% M = inv(P0)*P(t0+T) = Pf with
Pf=[solp(end,1),solp(end,2);solp(end,3),solp(end,4)];

% Monodromy matrix
M = Pf;
fprintf('Monodromy matrix:\n')
M
% Floquet multipliers
fprintf('Floquet multipliers:\n')
m=eig(Pf)
% Floquet exponents
B = 1/T*log(Pf);
fprintf('Floquet exponents:\n')
k = eig(B)
% Lyapunov exponents
fprintf('Lyapunov exponents:\n')
l = real(eig(B))
end

function rhs = vdp(t,x,p)
 rhs = [x(2);
        -p.mu*(x(1)^2-1)*x(2)-x(1)];
end

function rhs = fundeq(t,x,p)
 % Interpolation of the limit cycle at time t
 gam1 = interp1(p.tf,p.gam(:,1),t);
 gam2 = interp1(p.tf,p.gam(:,2),t);
 
 A = [0, 1; -2*p.mu*gam1*gam2-1, -p.mu*(gam1^2-1)];
 P = [x(1),x(2);x(3),x(4)];
 Pdot =A*P;
 rhs = [Pdot(1,1);Pdot(1,2);Pdot(2,1);Pdot(2,2)];
end	

sol_ex54

Ex5_4.m Objective-C source code, 2 KB

File contents

function FloquetVanDerPol
clc
fprintf('Determination of the Floquet and Lyapunov exponents for the stable\n')
fprintf('limit cycle of the van der Pol oscillator\n')
fprintf('       ddt x + mu*(x^2-1)*dt x + x = 0\n')
fprintf('with mu = 2.\n')
fprintf('\n')
% ------------------
% A. Schaum, SS-2017
% ------------------
close all
% Solution of the van der Pol oscillator dynamics
t = linspace(0,100,1000);
x0 = [0.1,0.2];
p.mu = 2;
[tv,solv] = ode45(@(t,x)vdp(t,x,p),t,x0);
figure(1)
plot(tv,solv(:,1),'LineWidth',2);

% Graphically determined starting and end time of one cylce
t0 = 84.48;%38.69;%31.17;%30.90455 
t1 = 92.19;%46.48;%38.69;%38.4425;
fprintf('Starting and end point of evaluation: t0 = %d, t1=%d\n\n', t0,t1)

%and period 
T  = t1-t0;
fprintf('Approximate period: T = %d\n\n', T)

% Interpolation of numerical solution over the interval [t0, t0+T]
p.tf = linspace(t0,t0+T,2000);
% Graphical comparison of the interpolation with the original solution
gamma = interp1(tv,solv,p.tf);
p.gam = gamma;
save gamma.mat gamma
figure(1);hold on;
plot(p.tf,p.gam(:,1),'r','LineWidth',2)

% Solution of the fundamental equation dt P=AP with A=df/dx evaluated at
% the approximated limit cycle p.gam.
% The matrix P = [p1,p2;p3,p4] is handled using the vector of its entries
% xp = [p1,p2,p3,p4]'. As starting point the identity matrix is chosen.
P0 = eye(2,2);
xp0 = [P0(1,1),P0(1,2),P0(2,1),P0(2,2)]';

[tp,solp] = ode45(@(t,x)fundeq(t,x,p),p.tf,xp0);

% The monodromy matrix is given by
% M = inv(P0)*P(t0+T) = Pf with
Pf=[solp(end,1),solp(end,2);solp(end,3),solp(end,4)];
M = Pf;
fprintf('Monodromy matrix:\n')
M
% Floquet multipliers
fprintf('Floquet multipliers:\n')
m=eig(Pf)
% Floquet exponents
B = 1/T*log(Pf);
fprintf('Floquet exponents:\n')
k = eig(B)
% Lyapunov exponents
fprintf('Lyapunov exponents:\n')
l = real(eig(B))
end

function rhs = vdp(t,x,p)
 rhs = [x(2);
        -p.mu*(x(1)^2-1)*x(2)-x(1)];
end

function rhs = fundeq(t,x,p)
 % Interpolation of the limit cycle at time t
 gam1 = interp1(p.tf,p.gam(:,1),t);
 gam2 = interp1(p.tf,p.gam(:,2),t);
 
 A = [0, 1; -2*p.mu*gam1*gam2 - 1, -p.mu*(gam1^2-1)];
 P = [x(1),x(2);x(3),x(4)];
 Pdot =A*P;
 rhs = [Pdot(1,1);Pdot(1,2);Pdot(2,1);Pdot(2,2)];
end	
Research