CoCalc Public Filesassignment4-octave-template (1).ipynbOpen with one click!
Author: Calum Henderson
Views : 63
Description: Jupyter notebook assignment4-octave-template (1).ipynb
Compute Environment: Ubuntu 18.04 (Deprecated)

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