The susceptible-infective-removed (SIR) epidemic model, where the population is
ID: 3111225 • Letter: T
Question
The susceptible-infective-removed (SIR) epidemic model, where the population is divided into susceptible S (t), infected I (t), and recovered R (t), in deterministic form is given by dS(t)/dt = -alpha I(t) S(t)/N dI(t)/dt = alpha I(s) S(t)/N - beta I(t) dR(t)/dt = beta I(t) where N = S (t) + I (t) + R(t) = S(0) + I(0) + R(0). The corresponding stochastic model is given by [4, 5, 67] [ds dI] = [- alpha IS/N alpha IS/N - beta I]dt + Squareroot alphaSI/N/Squareroot 2 + beta N/alpha S + 2 Squareroot beta N/alpha S [1 + beta N/alpha S - 1 -1 1 + beta N/alpha S + Squareroot beta N/alpha S][dW_1 dW_2] where W_1 and W_2 are the two Wiener processes. Taking alpha = 0.04, beta= 0.01, S (0) = 950, I (0) = 50 and time period [0, 100], obtain the numerical solution of the system and comment on the result.Explanation / Answer
In Matlab, the basic objects are matrices, i.e. arrays of numbers. Vectors can be thought of as special matrices. A row vector is recorded as a 1 × n matrix and a column vector is recorded as a m × 1 matrix. To enter a row vector in Matlab, type the following at the prompt ( > ) in the command window: > v = [0 1 2 3] and press enter. Matlab will print out the row vector. To enter a column vector type > u = [9; 10; 11; 12; 13] You can access an entry in a vector with > u(2) and change the value of that entry with > u(2)=47 You can extract a slice out of a vector with > u(2:4) You can change a column vector into a row vector, and row vector into a colomn vector easily in Matlab using > w = v’ (This is called transposing the vector and we call ’ the transpose operator.) There are also useful shortcuts to make vectors such as > x = -1:.1:1 and > y = linspace(0,1,11) Basic Formatting To make Matlab put fewer blank lines in its output, enter > format compact To make Matlab display more digits, enter > format long Note that this does not change the number of digits Matlab is using in its calculations; it only changes what is printed. Plotting Data Consider the data in Table 1.1.1 We can enter this data into Matlab with the following commands entered in the command window:
Table 1.1: Viscosity of a liquid as a function of temperature. > x = [ 5 20 30 50 55 ] > y = [ 0.08 0.015 0.009 0.006 0.0055] Entering the name of the variable retrieves its current values. For instance > x > y We can plot data in the form of vectors using the plot command: > plot(x,y) This will produce a graph with the data points connected by lines. If you would prefer that the data points be represented by symbols you can do so. For instance > plot(x,y,’*’) > plot(x,y,’o’) > plot(x,y,’.’) Data as a Representation of a Function A major theme in this course is that often we are interested in a certain function y = f(x), but the only information we have about this function is a discrete set of data {(xi , yi)}. Plotting the data, as we did above, can be thought of envisioning the function using just the data. We will find later that we can also do other things with the function, like differentiating and integrating, just using the available data. Numerical methods, the topic of this course, means doing mathematics by computer. Since a computer can only store a finite amount of information, we will almost always be working with a finite, discrete set of values of the function (data), rather than a formula for the function. Built-in Functions If we wish to deal with formulas for functions, Matlab contains a number of built-in functions, including all the usual functions, such as sin( ), exp( ), etc.. The meaning of most of these is clear. The dependent variable (input) always goes in parentheses in Matlab. For instance > sin(pi) should return the value of sin , which is of course 0 and > exp(0) will return e 0 which is 1. More importantly, the built-in functions can operate not only on single numbers but on vectors. For example > x = linspace(0,2*pi,40) > y = sin(x) > plot(x,y) will return a plot of sin x on the interval [0, 2] Some of the built-in functions in Matlab include: cos( ), tan( ), sinh( ), cosh( ), log( ) (natural logarithm), log10( ) (log base 10), asin( ) (inverse sine), acos( ), atan( ). To find out more about a function, use the help command; try > help plot
User-Defined Inline Functions If we wish to deal with a function that is a combination of the built-in functions, Matlab has a couple of ways for the user to define functions. One that we will use a lot is the inline function, which is a way to define a function in the command window. The following is a typical inline function: > f = inline(’2*x.^2 - 3*x + 1’,’x’) This produces the function f(x) = 2x 2 3x + 1. To obtain a single value of this function enter > f(2.23572) Just as for built-in functions, the function f as we defined it can operate not only on single numbers but on vectors. Try the following: > x = -2:.2:2 > y = f(x) This is an example of vectorization, i.e. putting several numbers into a vector and treating the vector all at once, rather than one component at a time, and is one of the strengths of Matlab. The reason f(x) works when x is a vector is because we represented x 2 by x.^2. The . turns the exponent operator ^ into entry-wise exponentiation, so that [-2 -1.8 -1.6].^2 means [(2)2 ,(1.8)2 ,(1.6)2 ] and yields [4 3.24 2.56]. In contrast, [-2 -1.8 -1.6]^2 means the matrix product [2, 1.8, 1.6][2, 1.8, 1.6] and yields only an error. The . is needed in .^, .*, and ./. It is not needed when you * or / by a scalar or for +. The results can be plotted using the plot command, just as for data: > plot(x,y) Notice that before plotting the function, we in effect converted it into data. Plotting on any machine always requires this step.
Function Programs Begin by clicking on the new document icon in the top left of the Matlab window (it looks like an empty sheet of paper). In the document window type the following: function y = myfunc ( x) y = 2* x .^2 - 3* x + 1; Save this file as: myfunc.m in your working directory. This file can now be used in the command window just like any predefined Matlab function; in the command window enter: > x = -2:.1:2; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Produces a vector of x values. > y = myfunc(x); . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Produces a vector of y values. > plot(x,y) Note that the fact we used x and y in both the function program and in the command window was just a coincidence. In fact, it is the name of the file myfunc.m that actually mattered, not what anything in it was called. We could just as well have made the function function nonsense = yourfunc ( inputvector ) nonsense = 2* inputvector .^2 - 3* inputvector + 1; Look back at the program. All function programs are like this one, the essential elements are: • Begin with the word function. • There is an input and an output. • The output, name of the function and the input must appear in the first line. • The body of the program must assign a value to the output. Functions can have multiple inputs, which are separated by commas. For example: function y = myfunc2d (x ,p ) y = 2* x .^ p - 3* x + 1; Functions can have multiple outputs, which are collected into a vector. Open a new document and type: function [ x2 x3 x4 ] = mypowers (x) x2 = x .^2; x3 = x .^3; x4 = x .^4; 5 6 LECTURE 2. MATLAB PROGRAMS Save this file as mypowers.m. In the command window, we can use the results of the program to make graphs: > x = -1:.1:1 > [x2 x3 x4] = mypowers(x); > plot(x,x,’black’,x,x2,’blue’,x,x3,’green’,x,x4,’red’) Script Programs Matlab uses a second type of program that differs from a function program in several ways, namely: • There are no inputs and outputs. • A script program may use and change variables in the current workspace (the variables used by the command window.) Below is a script program that accomplishes the same thing as the function program plus the commands in the previous section: x2 = x .^2; x3 = x .^3; x4 = x .^4; plot (x ,x , ’ black ’ ,x ,x2 , ’ blue ’,x ,x3 , ’ green ’,x ,x4 , ’ red ’) Type this program into a new document and save it as mygraphs.m. In the command window enter: > x = -1:.1:1; > mygraphs Note that the program used the variable x in its calculations, even though x was defined in the command window, not in the program. Many people use script programs for routine calculations that would require typing more than one command in the command window. They do this because correcting mistakes is easier in a program than in the command window. Program Comments For programs that have more than a couple of lines it is important to include comments. Comments allow other people to know what your program does and they also remind yourself what your program does if you set it aside and come back to it later. It is best to include comments not only at the top of a program, but also with each section. In Matlab anything that comes in a line after a % is a comment. For a function program, the comments should at least give the purpose, inputs, and outputs. A properly commented version of the function with which we started this section is: function y = myfunc ( x) % Computes the function 2 x ^2 -3 x +1 % Input : x -- a number or vector ; % for a vector the computation is elementwise % Output : y -- a number or vector of the same size as x y = 2* x .^2 - 3* x + 1;
Related Questions
drjack9650@gmail.com
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.