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	

ex61.py

ex6_1.py Python Source, 2 KB

File contents

#------------------------------------------------------------
# Solution file for the numerical simulation of exercise 6.1 
# Applied Nonlinear Dynamics  etit-614, SS2017
# A. Schaum, 02.06.2017
#------------------------------------------------------------

import numpy as np
from scipy.integrate import odeint, ode
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation

def ddt(x,t):
    	# right-hand side of ode
	x1 = x[0]
	x2 = x[1]
	dx1dt = x1-x2-x1*(x1**2+x2**2)
	dx2dt = x1+x2-x2*(x1**2+x2**2)
	return [dx1dt, dx2dt]

def make_line(x1,x2):
 	# creating 3d line from simulation data
	V = x1**2+x2**2
	line = np.vstack((x1,x2,V))
	return line

def update_lines(num, dataLines, lines):
    for line, data in zip(lines, dataLines):
        line.set_data(data[0:2, num-1:num])
        line.set_3d_properties(data[2,num-1:num])
    return lines

# Solving the differential equation
# 1. Creating a time vector
dt = 0.05
t = np.arange(0.0, 50, dt)
# 2. Setting the initial condition
x0 = np.array([-2,2])
# 4. Solving the ode with odeint
sol = odeint(ddt,x0,t)

x1=sol[:,0]
x2=sol[:,1]

# Creating the line object
data = make_line(x1,x2)

# Creating 3d plot object
fig = plt.figure()
ax = p3.Axes3D(fig)

# Setting the axes properties
ax.set_xlim3d([-2.0, 2.0])
ax.set_xlabel('x_1')

ax.set_ylim3d([-2.0, 2.0])
ax.set_ylabel('x_2')

ax.set_zlim3d([0.0, 20.0])
ax.set_zlabel('V')

ax.set_title('Simulation for exercise 6.1 (AND-SS2017)')

x1pl = np.linspace(-2.0,2.0,100)
x2pl = np.linspace(-2.0,2.0,100)

# Creating meshgrid and energy well
X, Y = np.meshgrid(x1pl,x2pl)

V = X**2+Y**2
# Plotting the energy well
p = ax.plot_surface(X,Y,V,rstride=2,cstride=2)
# Setting the view angle
ax.view_init(elev=50., azim=30.)
# Creating plot data
data = [make_line(x1,x2)]

lines = [ax.plot(data[0][0,0:1], data[0][1,0:1], data[0][2,0:1], 'o')[0]]
# Animating the simulation result
ani = animation.FuncAnimation(fig, update_lines, np.arange(1, len(x1)),
                              fargs=(data, lines), interval=20, blit=True)
# Saving the animation as mp4 (takes some time..can be commented out for demonstration or test purposes

#ani.save('ex6_1.mp4', fps=20)

# Show simulation results
plt.show()

Research