CEIC3000 Assignment 4¶

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

pkg load odepkg

Question 1 template¶

Part (a)¶

% 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)¶

% 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 %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 %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)¶

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)¶

figure
hold on

% define a time range
time = linspace(0,20,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('\theta','FontSize',7)
ylabel('d\theta/dt','FontSize',7)

title(['Phase Portrait for the Damping of Pendulous Motion'...
,' at Angle \theta'],'FontSize',6)
axis('square')

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

hold off

Question 2 template¶

Part (a)¶

% 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)

% 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

% define a time range
time = linspace(0,10,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)¶

% 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,10,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