ex1sol

exercise1_sol.pdf PDF document, 82 KB

ex1_2_code

solveDiffEq.py Python Source, 1 KB

File contents

#-----------------------------------
# Applied Nonlinear dynamics
# Solution of exercise 1.2
#-----------------------------------
# Python code for solving a nonlinear continuous-time differential equation
#-----------------------------------

# Import relevant python libraries
import numpy as np
import scipy.linalg as linalg
from scipy.integrate import ode
import matplotlib.pyplot as plt

# Definition of the rhs-function
def dxdt(t,x,a):
  return a*x*(1.0-x**2.0) #returns the right-hand side of the diff. eq.

#--------------------------------------
# Preparation of the necessary variables containing system structure and parameters, as well as simulation parameters
#--------------------------------------

# Simulation parameters:
# Initial condition
x0 = 0.5
# System parameter
a = 2.0
# Simulation time
Tsim = 10.0

#----------------------------------------
# Solution of the continuous-time model
#----------------------------------------
# Setting the structure of the simulation object according tp scipy.integrate.ode standard
ss = ode(dxdt, jac=None) #No explicit Jacobian of the dynamics used
ss.set_integrator('lsoda',method='bdf',max_step=0.001)#alternatively 'dopri5' for 5th order Runge-Kutta method
ss.set_initial_value(x0)
ss.set_f_params(a)

# Determining the solution between prescribed points (the ones later used for the graphical output)
i = 1
ts = np.linspace(0,Tsim,500)
dt = ts[1]-ts[0]
# Declaring the state variable
xs = np.zeros(shape = (1,ts.size))
# Setting the initial condition
xs[:,0] = x0
# Loop for solution over all time intervals in ts-vector
while ss.successful() and i < ts.size:
    ss.integrate(ss.t + dt)
    xs[:,i] = ss.y
    i += 1 

#--------------------------------------
# Graphical output
#--------------------------------------
fig = plt.figure(figsize=(6,6))
ax1 = fig.add_subplot(1,1,1)

p = ax1.plot(ts,xs[0,:])
ax1.grid()
ax1.set_ylabel(r'$x_1$')
ax1.set_xlabel(r'$t$')

plt.show() #Needed to show the figure..

ex33.m

Ex33sol.m Objective-C source code, 2 KB

File contents

function Ex33sol(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 Ex33sol(x0)
	% where x0 is the initial condition for the simulation.
	% A. Schaum, 2019


	% Simulation time vector
	T = 2*pi % 1 period of the matrix A
	t = linspace(0,3*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 = [-0.5+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

Research