CEIC3000 Assignment 4

In [1]:
% Name: [Calum Henderson]
% zID: [z3464226]
Out[1]:
In [2]:
pkg load odepkg
Out[2]:

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])
Out[3]:
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 %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)))
Out[4]:
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
Gnuplot Produced by GNUPLOT 4.6 patchlevel 6 -2 0 2 4 6 8 10 -5 0 5 10 15 x2 x1 Trajectories Starting from x1 = p/16 and x2 = 1, 4, 9 gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a
Out[5]:

Part (e)

In [6]:
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
Gnuplot Produced by GNUPLOT 4.6 patchlevel 6 -4 -3 -2 -1 0 1 2 3 4 -2 0 2 4 6 8 10 12 dq/dt q Phase Portrait for the Damping of Pendulous Motion at Angle q gnuplot_plot_1a gnuplot_plot_2a gnuplot_plot_3a gnuplot_plot_4a gnuplot_plot_5a gnuplot_plot_6a gnuplot_plot_7a gnuplot_plot_8a gnuplot_plot_9a gnuplot_plot_10a gnuplot_plot_11a gnuplot_plot_12a gnuplot_plot_13a gnuplot_plot_14a gnuplot_plot_15a gnuplot_plot_16a gnuplot_plot_17a gnuplot_plot_18a gnuplot_plot_19a gnuplot_plot_20a gnuplot_plot_21a gnuplot_plot_22a gnuplot_plot_23a gnuplot_plot_24a gnuplot_plot_25a gnuplot_plot_26a gnuplot_plot_27a gnuplot_plot_28a gnuplot_plot_29a gnuplot_plot_30a gnuplot_plot_31a gnuplot_plot_32a gnuplot_plot_33a gnuplot_plot_34a gnuplot_plot_35a gnuplot_plot_36a gnuplot_plot_37a gnuplot_plot_38a gnuplot_plot_39a gnuplot_plot_40a gnuplot_plot_41a gnuplot_plot_42a gnuplot_plot_43a gnuplot_plot_44a gnuplot_plot_45a gnuplot_plot_46a gnuplot_plot_47a gnuplot_plot_48a gnuplot_plot_49a gnuplot_plot_50a gnuplot_plot_51a gnuplot_plot_52a gnuplot_plot_53a gnuplot_plot_54a gnuplot_plot_55a gnuplot_plot_56a gnuplot_plot_57a gnuplot_plot_58a gnuplot_plot_59a gnuplot_plot_60a gnuplot_plot_61a gnuplot_plot_62a gnuplot_plot_63a gnuplot_plot_64a gnuplot_plot_65a gnuplot_plot_66a gnuplot_plot_67a gnuplot_plot_68a gnuplot_plot_69a gnuplot_plot_70a gnuplot_plot_71a gnuplot_plot_72a gnuplot_plot_73a gnuplot_plot_74a gnuplot_plot_75a gnuplot_plot_76a gnuplot_plot_77a gnuplot_plot_78a gnuplot_plot_79a gnuplot_plot_80a gnuplot_plot_81a gnuplot_plot_82a gnuplot_plot_83a gnuplot_plot_84a gnuplot_plot_85a gnuplot_plot_86a gnuplot_plot_87a gnuplot_plot_88a gnuplot_plot_89a gnuplot_plot_90a gnuplot_plot_91a gnuplot_plot_92a gnuplot_plot_93a gnuplot_plot_94a gnuplot_plot_95a gnuplot_plot_96a gnuplot_plot_97a gnuplot_plot_98a gnuplot_plot_99a gnuplot_plot_100a
Out[6]:

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
Out[7]:
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 [ ]:
% 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)))
Out[ ]:
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 [ ]:
% 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)

In [ ]:
% 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