CoCalc Shared Filesassignment4-octave-template (1).ipynb
Author: Calum Henderson
Views : 4
Description: Jupyter notebook assignment4-octave-template (1).ipynb

# CEIC3000 Assignment 4

In [1]:
% Name: [Calum Henderson]
% zID: [z3464226]

In [2]:
pkg load odepkg


## Question 1 template

### Part (a)

In [3]:
% define parameters
k = 0.8;

% define the model
dxdt = @(t,x) [x(2); -sin(x(1))-k*x(2)];

% confirm steady state 1 is at (0,0)
xe1 = fsolve(@(x) dxdt(0,x),[0,0])

% confirm steady state 2 is at (pi,0)
xe2 = fsolve(@(x) dxdt(0,x),[pi,0])

xe1 = 0 0 xe2 = 3.14159 0.00000

### Part (b)

In [4]:
% define the Jacobian
Jac = @(xe) [0 1; -cos(xe(1)) -k];

% linearise the model around steady state 1
A1 = Jac(xe1);

% linearise the model around steady state 2
A2 = Jac(xe2);

% determine the stability of steady state 1 and 2
[V1 D1] = eig(A1);
[V2 D2] = eig(A2);

% print the eigenvalues of the model linearised around (0,0)
fprintf(['eigenvalues of the Jacobian evaluated'...
,' at the (0,0) steady state are \n %s and %s \n \n']...
, num2str(D1(1,1)), num2str(D1(2,2)))

% print the eigenvalues of the model linearised around (pi,0)
fprintf(['eigenvalues of the Jacobian evaluated'...
,' at the (pi,0) steady state are \n %s and %s \n \n']...
, num2str(D2(1,1)), num2str(D2(2,2)))

eigenvalues of the Jacobian evaluated at the (0,0) steady state are -0.4+0.91652i and -0.4-0.91652i eigenvalues of the Jacobian evaluated at the (pi,0) steady state are -1.477 and 0.67703

### Part (d)

In [5]:
figure
hold on

% define a time range
time = linspace(0,50,1001);

% solve the ODE from (-pi/16,1)
[t,x] = ode45(dxdt,time,[-pi/16;1]);

% plot the trajectory starting at (-pi/16,1)
plot(x(:,1),x(:,2),'r')
hold on

% solve the ODE from (-pi/16,4)
[t,x] = ode45(dxdt,time,[-pi/16,4]);

% plot the trajectory starting at (-pi/16,4)
plot(x(:,1),x(:,2),'b')
hold on

% solve the ODE from (-pi/16,9)
[t,x] = ode45(dxdt,time,[-pi/16,9]);

% plot the trajectory starting at (-pi/16,9)
plot(x(:,1),x(:,2),'g')

% make the plot look good
xlabel('x_{1}');
ylabel('x_{2}');

title(['Trajectories Starting from x_{1} = pi/16 and'...
,' x_{2} = 1, 4, 9'],'FontSize',7);

axis('square');

xt = get(gca, 'XTick');
set(gca, 'FontSize', 6);

hold off


### Part (e)

In [6]:
figure
hold on

% define a time range
time = linspace(0,100,1001);

% define the the number of trajectories to be plotted
n_xs = 100;

% define the starting positions of the trajectories
% from the bottom of the ROI
x1s = linspace(-2*pi,6*pi,n_xs/2);
x2s = ones(1,n_xs/2).*-4;

% create a matrix of starting points for the trajectories
% along the bottom of the ROI
x0s = [x1s; x2s];

% iterate over each x0 and draw in the trajectory
% (make sure x0s has the right shape for the iteration)
for x0 = x0s
[t,x] = ode45(dxdt,time,x0);
plot(x(:,1),x(:,2))
hold on
end

% redefine the starting positions of the trajectories
% to start from the top of the ROI
x2s = ones(1,n_xs/2).*4;

% recreate the matrix of starting points for the trajectories
% along the top of the ROI
x0s = [x1s; x2s];

% iterate over each x0 and draw in the trajectory
% (make sure x0s has the right shape for the iteration)
for x0 = x0s
[t,x] = ode45(dxdt,time,x0);
plot(x(:,1),x(:,2))
hold on
end

% make plot look good
xlim([-pi, 4*pi])
ylim([-4, 4])
xlabel('x_{1}','FontSize',7)
ylabel('dx_{1}/dt','FontSize',7)

title(['Phase Portrait for the Damping of Pendulous Motion'...
,' at Angle x_{1}'],'FontSize',6)
axis('square')

xt = get(gca, 'XTick');
set(gca, 'FontSize', 6);

hold off


## Question 2 template

### Part (a)

In [7]:
% define parameters
F = 0.5;
x1i = 1.90;
k = 0.1;

% define the model
dxdt = @(t,x) [F*(x1i - x(1))-x(1)*x(2)^2;
-(k+F)*x(2)+x(1)*x(2)^2];

% solve for each x2 steady state by
% determining the roots of the cubic polnomial:
% -1.2*x^3 + 1.9*x^2 - 0.6*x = 0
p = [-1.2, 1.9, -0.6, 0];
x2 = roots(p);

% solve for each x1 steady state by
% substituting the x2 values into
% the formula found analytically
x1 = 1.9 - 1.2.*x2;

% confirm each steady state with the fsolve function
xe1 = fsolve(@(x) dxdt(0,x),[x1(1),x2(1)]);
xe2 = fsolve(@(x) dxdt(0,x),[x1(2),x2(2)]);
xe3 = fsolve(@(x) dxdt(0,x),[x1(3),x2(3)]);

% print the steady state points
if [x1(1), x2(1)] == xe1
fprintf(['steady state 1 is located at (%s, %s) \n \n']...
, num2str(x1(1)), num2str(x2(1)))
else fprintf(['there is not a steady state at (%s, %s) \n \n']...
, num2str(x1(1)),num2str(x2(1)))
end

if [x1(2), x2(2)] == xe2
fprintf(['steady state 2 is located at (%s, %s) \n \n']...
, num2str(x1(2)), num2str(x2(2)))
else fprintf(['there is not a steady state at (%s, %s) \n \n']...
, num2str(x1(2)),num2str(x2(2)))
end

if [x1(3), x2(3)] == xe3
fprintf(['steady state 3 is located at (%s, %s) \n \n']...
, num2str(x1(3)), num2str(x2(3)))
else fprintf(['there is not a steady state at (%s, %s) \n \n']...
, num2str(x1(3)),num2str(x2(3)))
end

steady state 1 is located at (0.5228, 1.1477) steady state 2 is located at (1.3772, 0.43567) steady state 3 is located at (1.9, 0)
In [8]:
% define the Jacobian
Jac = @(xe) [-F-xe(2)^2, -2*xe(1)*xe(2);
xe(2)^2, -(k+F)+2*xe(1)*xe(2)];

% linearise the model around steady state 1
A1 = Jac(xe1);

% linearise the model around steady state 2
A2 = Jac(xe2);

% linearise the model around steady state 2
A3 = Jac(xe3);

% determine the stability of steady state 1 and 2
[V1 D1] = eig(A1);
[V2 D2] = eig(A2);
[V3 D3] = eig(A3);

% print the eigenvalues
fprintf(['eigenvalues of the Jacobian evaluated at'...
,' steady state 1 are \n %s and %s \n \n']
, num2str(D1(1,1)), num2str(D1(2,2)))
fprintf(['eigenvalues of the Jacobian evaluated at'...
,' steady state 2 are \n %s and %s \n \n']
, num2str(D2(1,1)), num2str(D2(2,2)))
fprintf(['eigenvalues of the Jacobian evaluated at'...
,' steady state 3 are \n %s and %s \n \n']
, num2str(D3(1,1)), num2str(D3(2,2)))

eigenvalues of the Jacobian evaluated at steady state 1 are -0.60857+0.3463i and -0.60857-0.3463i eigenvalues of the Jacobian evaluated at steady state 2 are -0.47865 and 0.38884 eigenvalues of the Jacobian evaluated at steady state 3 are -0.6 and -0.5

### Part (b)

Hint: draw trajectories starting on the left and right sides of the figure

In [9]:
% define a time range
time = linspace(0,100,1001);

% define the number of trajectories to be plotted
n_xs = 100;

% define the starting positions of the trajectories
% from the left side of the ROI
x2s = linspace(0,3,n_xs/2);
x1s = zeros(1,n_xs/2);

% create a matrix of starting points for the trajectories
% along the left side of the ROI
x0s = [x1s;x2s];

% iterate over each x0 and draw in the trajectory
% (make sure x0s has the right shape for the iteration)
for x0 = x0s
[t,x] = ode45(dxdt,time,x0);
plot(x(:,1),x(:,2))
hold on
end

% redefine the starting positions of the trajectories
% to start from the right side of the ROI
x1s = ones(1,n_xs/2).*3;

% recreate the matrix of starting points for the trajectories
% along the right side of the ROI
x0s = [x1s;x2s];

% iterate over each x0 and draw in the trajectory
% (make sure x0s has the right shape for the iteration)
for x0 = x0s
[t,x] = ode45(dxdt,time,x0);
plot(x(:,1),x(:,2))
hold on
end

% make the plot look good
xlim([0, 3])
ylim([0, 3])
xlabel('C_{A}','FontSize', 7)
ylabel('C_{B}','FontSize', 7)
title(['Phase Portrait of the Autocatalytic Reaction'...
,' Occuring in a CSTR'], 'FontSize', 6)

xt = get(gca, 'XTick');
set(gca, 'FontSize', 6);
axis('square')

hold off


### Part (c)

In [10]:
% define the new parameter
x1i = 1.50;

% redefine the model with the new parameter
dxdt = @(t,x) [F*(x1i - x(1))-x(1)*x(2)^2;
-(k+F)*x(2)+x(1)*x(2)^2];

% define a time range
time = linspace(0,100,1001);

% define the number of trajectories to be plotted
n_xs = 100;

% define the starting positions of the trajectories
% from the left side of the ROI
x2s = linspace(0,3,n_xs/2);
x1s = zeros(1,n_xs/2);

% create a matrix of starting points for the trajectories
% along the left side of the ROI
x0s = [x1s;x2s];

% iterate over each x0 and draw in the trajectory
% (make sure x0s has the right shape for the iteration)
for x0 = x0s
[t,x] = ode45(dxdt,time,x0);
plot(x(:,1),x(:,2))
hold on
end

% redefine the starting positions of the trajectories
% to start from the right side of the ROI
x1s = ones(1,n_xs/2).*3;

% recreate the matrix of starting points for the trajectories
% along the right side of the ROI
x0s = [x1s;x2s];

% iterate over each x0 and draw in the trajectory
% (make sure x0s has the right shape for the iteration)
for x0 = x0s
[t,x] = ode45(dxdt,time,x0);
plot(x(:,1),x(:,2))
hold on
end

% make the plot look good
xlim([0, 3])
ylim([0, 3])
xlabel('C_{A}','FontSize', 7)
ylabel('C_{B}','FontSize', 7)
title(['Phase Portrait of the Autocatalytic Reaction'...
,' Occuring in a CSTR with x_{1i} = 1.50'], 'FontSize', 5.1)

xt = get(gca, 'XTick');
set(gca, 'FontSize', 6);
axis('square')

hold off