Differential Equations An Introduction for Engineers

2 + 3

ans = 5

2 + 3;

>> 4*6 
ans = 24

[1 2 3] * [4 3 2];

>> [1 2 3] .* [4 3 2]
ans = 
   [4 6 6]
>> [1 4 6].^2
ans = 
    [1 16 36]
>> [5 4 2] ./ [10 2 6]
ans = 
    [0.5 2 0.3333]

val = 2+3;

>> val * 4
ans = 
    20
    
>> val^2 + 2
ans = 
    27

ans = 5

numVar = sqrt(15);
strVar = ``Hello World!'';

A = [2,4,6];
B = [1 3 5];

C = [5;7;8];
M = [1,2,3;5,6,7];

>> M = [1,2,3;5,6,7];
>> M(1,1)
ans =
     1

>> M(1,3)
ans =
     3

>> M(2,1)
ans =
     5

>> M(3,1) 
Index in position 1 exceeds array bounds (must not exceed 2).

function [a,b] = testFunction(x, y, z)
    % Code here
end

function [a,b] = testFunction(x, y, z)
    a = x+y+z;
    b = x*y*z;
end

function [a,b] = testFunction(x, y, z)
    a = x+y+z;
    b = x.*y.*z;
end

f = @(x,y) x.^2 + 4.*x.*y + y.^2;

>> f(2,1)
ans = 
    13

fShort = @(x,y) x.^2 + y.^2;

function z = fLong(x,y)
z = x.^2 + y.^2;
end

ode45(fShort, [0, 3], 1)

ode45(fLong, [0, 3], 1)

ode45(@fLong, [0, 3], 1)

fNew = @(y) fShort(1,y)

xPts = [1,2,3,4,5];
fx = @(x) x.^2 + 2;
yPts = [2,3,2,3,1];
figure(1);
plot(xPts, yPts);
figure(2);
plot(xPts, fx(xPts));

xPts = linspace(1,5,100);
fx = @(x) x.^2 + 2;
gx = @(x) x.^2 - 3*x + 7;
figure(1);
hold on;
plot(xPts, fx(xPts));
plot(xPts, gx(xPts));
hold off;

function quiver244(f, t_min, t_max, y_min, y_max, col)
% quiver244.m
% Author: Matt Charnley
%
% This function draws a quiver plot for the ODE dy/dt = f(t,y) for
% t_min <= t <= t_max and y_min <= y <= y_max. The function f should be
% passed in as an anonymous function, of two variables or as a function
% handle
%
% The function draws this quiver plot in color col and saves it on the 
% current figure, and generates a normalized version
% (all vectors are the same length) as the next figure,
% so that it can be accessed outside of this function.
% For this second figure, the magnitude of the arrows does not mean
% anything, but it is easier to see the direction of them. 
% so that it can be accessed outside of this function. It will start with
% hold on; and end with hold off;, so the figure needs to be cleared in the
% main file if needed.

f = @(t,y) t - exp(y);
quiver244(f, 0, 5, -6, 6, 'b');

quiver244(@f2, 0, 5, -6, 6, 'b');

function z = f2(t,y)
    z = t - exp(y);
end

function samplePlots244(f, t_min, t_max, y_min, y_max, t_0, y_0, col)
% This function takes the ODE dy/dt = f(t,y) and plots sample solutions
% with initial value (t_0, y_0). It uses ode45 to sketch out the solutions.
% t_0 must be between t_min and t_max. It also truncates the function f so
% that functions will not go off to infinity, causing this to work properly
% on vector inputs for initial conditions in y. The input y_0 can be a vector 
% of initial values, and this function will plot a curve
% for each of those values. If using a vector of initial
% conditions, the function must be written with vector element-wise
% operations.

f = @(t,y) y.*(y-5).*(y+6);
samplePlots244(f, -1, 6, -7, 6, 0, [-1,0.5,4,5], 'r')

function bifDiag244(f, a_min, a_max, y_min, y_max)
% This function draws a bifurcation diagram for the ode dy/dt = f(alpha, y)
% with parameter alpha running from a_min to a_max. The axes are
% constrained to be from a_min to a_max in the horizontal direction and
% y_min to y_max in the vertical direction.
%
% The black marks are for equilibrium solutions, the blue regions are where
% the solution will tend upwards, and the red region is where it will tend
% downwards.

f = @(a,y) y.^2 - a.^2;
bifDiag244(f, -3, 3, -5, 5, 3);

function quiver2D244(f,g, x_min, x_max, y_min, y_max, col)
% quiver2D244.m
% Author: Matt Charnley
%
% This function draws a quiver plot for the ODE dx/dt = f(x,y), dy/dt = g(x,y) for
% x_min <= x <= x_max and y_min <= y <= y_max. The functions f and g should be
% passed in as an anonymous functions, f = @(x,y) ... 
%
% The function draws this quiver plot in color col in the current figure
% and generates a normalized version (all vectors are the same length)
% as the next figure, so that it can be accessed outside of this function.
% For this second figure, the magnitude of the arrows does not mean
% anything, but it is easier to see the direction of them. 
%
% It will start with
% hold on; and end with hold off;, so the figure needs to be cleared in the
% main file if needed.

f = @(x,y) 3.*x - 2.*x.*y;
g = @(x,y) 2.*y - 3.*x.*y;
quiver2D244(f,g, 0, 5, 0, 5, 'g');

function phaseLine(f, ymin, ymax)
% This function draws a representation of the phaseline for the
% differential equation dy/dt = f(y). The graph is drawn from ymin to ymax,
% and looks for solutions to f(y) = 0 in that region to find equilibrium
% solutions. This requires the Optimization Toolbox fsolve to run
% correctly.

f = @(y) y.*(y-3).*(y+2);
phaseLine(f, -4, 5);

function phasePortrait244(F, G, xmin, xmax, ymin, ymax, tmin, tmax, x0, y0)
% This function draws a 2 dimensional phase portrait for the system dx/dt =
% F(x,y) and dy/dt = G(x,y). The phase portrait will be draw with x bounds
% xmin <= x <= xmax and ymin <= y <= ymax. It is assumed that the initial 
% conditions x0 and y0 are at $t=0$, with tmin <= 0 and tmax >=0. x0 and y0
% can be inputted as vectors that are the same length, and a sample curve
% will be drawn for each of them. The black dot will always be plotted at tmin.

f = @(x,y) 2.*x - 3.* y;
g = @(x,y) -3.*x + y;
phasePortrait244(f, g, -3, 3, -3, 3, -2, 2, [1, 0, -1, 1, 0, -1], [1,1,1,-1,-1,-1]);

function [t, y] = rungeKuttaMethod(f, dt, Tf, T0, y0)
% This method solves the ODE dy/dt = f(t, y) using the Runge Kutta method
% from t=T0 to t = Tf with time step dt and initial condition y0 at t = T0.
% In this case, f should be a function of two variables, t
% (time) and y.

function [t,y] = rungeKuttaSystemMethod(f, T0, Tf, dt, y0)
% This method solves the ODE system dy/dt = f(t, y) using the Runge Kutta method
% from t=T0 to t = Tf with time step dt and initial condition y0 at t = T0.
% In this case, f should be a vector valued function of two variables, t
% (time) and y (n-dimensional vector of unknowns). The length of the vector
% y0 will determine the size of the system.

function [S,I,R] = SIRModel_244(r, c, ICs, Tf)
% This code runs an SIR model for disease spread. The system of differential equations used here is
%   S' = -r*S*I
%   I' = r*S*I - cI
%   R' = c*I
%
% The solution is computed using the RungeKutta method, with the helper
% method rungeKuttaSystemMethod. The system is solved from t=0 to t=Tf,
% with initial conditions ICs given as a 3 component vector.

function [S,I,Q,R,D] = SIRQModel_244(alpha, beta, gamma, delta, eta, rho, ICs, Tf)
% This code runs a more complicated SIR model that adds in Q (a quarantined
% population) and D (a deceased population). The system of differential equations used here is
%   S' = -alpha*S*I
%   I' = alpha*S*I - (beta+gamma+delta)I
%   Q' = beta*I - (eta + rho)Q
%   R' = gamma*I + eta*Q
%   D' = delta*I + rho*Q
%
% The solution is computed using the RungeKutta method, with the helper
% method rungeKuttaSystemMethod. The system is solved from t=0 to t=Tf,
% with initial conditions ICs given as a 5 component vector.

function [S,I,Q,R,D] = SIRQVModel_244(alpha, beta, gamma, delta, eta, rho, zeta, ICs, Tf)
% This code runs a more complicated SIR model that adds in Q (a quarantined
% population) and D (a deceased population). The V component adds
% vaccination into the picture, where members are moved from S to R
% directly. The system of differential equations used here is
%   S' = -alpha*S*I - zeta*S
%   I' = alpha*S*I - (beta+gamma+delta)I
%   Q' = beta*I - (eta + rho)Q
%   R' = gamma*I + eta*Q+zeta*S
%   D' = delta*I + rho*Q
%
% The solution is computed using the RungeKutta method, with the helper
% method rungeKuttaSystemMethod. The system is solved from t=0 to t=Tf,
% with initial conditions ICs given as a 5 component vector.

[S,I,R] = SIRModel_244(0.1, 0.2, [0.99; 0.01; 0], 400);
[S,I,Q,R,D] = SIRQModel_244(0.15, 0.08, 0.02, 0.03, 0.01, 0.04, [0.95; 0.05; 0; 0; 0], 400);
[S,I,Q,R,D] = SIRQVModel_244(0.15, 0.08, 0.02, 0.03, 0.01, 0.04,0.2, [0.95; 0.05; 0; 0; 0], 400);