Labs System Modeling

Lab 1: MATLAB DEVELOPMENT ENVIRONMENT

Objectives

Learn about Matlab Development Environment
Learn about the starting and quitting Matlab
Learn about Matlab desktop
Learn about Matlab desktop tools

Prerequisites

Students are expected to be familiar with Windows system
Students should have basic knowledge of any text editors

Requirement

Matlab software any version

Remarks

This lab session will be evaluated and it will carry 2 % of your activity credit points.
All the source codes and executables related to this lab session should be saved in a folder
named Lab1. If not, zero marks will be given. Before you leave the lab make sure you show to your instructor that you have saved your file in Lab1 folder.

Starting and Quitting MATLAB

Starting MATLAB

On Windows platforms, to start MATLAB, double-click the MATLAB shortcut icon on your Windows desktop.

You can change the directory in which MATLAB starts, define startup options including running a
script upon startup, and reduce startup time in some situations.

Quitting MATLAB

To end your MATLAB session, select Exit MATLAB from the File menu in the desktop, or type
quit” in the Command Window. To execute specified functions each time MATLAB quits, such as saving the workspace, you can create and run a name.m” script.

MATLAB Desktop

When you start MATLAB, the MATLAB desktop appears, containing tools (graphical user interfaces) for managing files, variables, and applications associated with MATLAB.

You can change the way your desktop looks by opening, closing, moving, and resizing the tools in it.
Use the View menu to open or close the tools. You can also move tools outside the desktop or move
them back into the desktop (docking). All the desktop tools provide common features such as context
menus and keyboard shortcuts.

You can specify certain characteristics for the desktop tools by selecting Preferences from the File
menu. For example, you can specify the font characteristics for Command Window text. For more
information, click the Help button in the Preferences dialog box.

Desktop Tools

This section provides an introduction to the MATLAB desktop tools. You can also use MATLAB
functions to perform most of the features found in the desktop tools. The tools are

  • Command Window
  • Command History
  • Start Button and Launch Pad
  • Help Browser
  • Current Directory Browser
  • Workspace Browser
  • Array Editor
  • Editor/Debugger

Command Window

Use the Command Window to enter variables and run functions and M-files. In
this example, MAGIC(4) is called to create an 4-by-4 matrix with the integers from 1 through 4 = 16
with equal row, column, and diagonal sums.

A=magic(4);

Command History

Statements you enter in the Command Window are logged in the Command History. In the Command
History, you can view previously run statements, and copy and execute selected statements. In green is given the date/time when the execution of the commands took place.

To save the input and output from a MATLAB session to a file, use the diary function as the follows:

diary filename

After finishing writing into the file, use

diary off

to close the file.

Start Button

The MATLAB Start button provides easy access to tools, demos, and documentation. Just click the
button to see the options as shown in Figure below.

fig1.png

Help Browser

Use the Help browser to search and view documentation and demos for all your MathWorks products.
The Help browser is a Web browser integrated into the MATLAB desktop that displays HTML
documents.

To open the Help browser, click the help button in the toolbar,

fig2.png

Press F1, or type helpbrowser in
the Command Window and the following image will appear (see Figure below).

fig3.png

The Help browser consists of two panes, the Help Navigator, which you use to find information, and
the display pane, where you view the information. Tabs in the Help Navigator pane provide different
ways to find documentation and demos. Drag the separator bar to adjust the width of the panes. View
documentation in the display pane. Use the close box to hide the pane.

Help Navigator

Use the Help Navigator to find information. It includes

  • Product filter—set the filter to show documentation only for the products you specify.
  • Contents tab—view the titles and tables of contents of documentation for your products.
  • Index tab—find specific index entries (selected keywords) in the MathWorks documentation for your products.
  • Demos tab—view and run demonstrations for your MathWorks products.
  • Search tab—look for a specific word or phrase in the documentation. To get help for a specific function, set the Search type to Function Name.
  • Favorites tab—view a list of links to documents you previously designated as favourites.

Display Pane

After finding documentation using the Help Navigator, view it in the display pane. While viewing the documentation, you can

  • Browse to other pages—Use the arrows at the tops and bottoms of the pages to move through the document, or use the back and forward buttons in the toolbar to go to previously viewed pages.
  • Bookmark pages—Click the Add to Favorites button in the toolbar.
  • Print pages—Click the print button in the toolbar.
  • Find a term in the page—Type a term in the Find in page field in the toolbar and click Go.

Other features available in the display pane are copying information, evaluating a selection, and viewing Web pages.

For More Help

In addition to the Help browser, you can use help functions. To get help for a specific function, use doc. For example, doc format displays documentation for the format function in the Help browser.

fig4.png

If you type help followed by the function name, for instance ``help diary”, a briefer form of the documentation appears in the Command Window (see Figure below).

Other means for getting help include contacting Technical Support (http://www.mathworks.com/support) and participating in the newsgroup for MATLAB users, comp.soft-sys.matlab.

Current Directory Browser

MATLAB file operations use the current directory and the search path as reference points. Any file you want to run must either be in the current directory or on the search path.

A quick way to view or change the current directory is by using the Current Directory field in the desktop toolbar as shown in figure.

fig5.png

To search for, view, open, and make changes to MATLAB-related directories and files, use the MATLAB Current Directory browser. Alternatively, you can use the functions

dir
cd

and

delete

Search Path

MATLAB uses a search path to find M-files and other MATLAB-related files, which are organized in directories on your file system. Any file you want to run in MATLAB must reside in the current directory or in a directory that is on the search path. Add the directories containing files you create to the MATLAB search path. By default, the files supplied with MATLAB and MathWorks toolboxes are included in the search path.
To see which directories are on the search path or to change the search path, select Search Path from the File menu in the desktop and use the Set Path dialog box (shown in the figure below). Alternatively, you can use the path function to view the search path, addpath to add directories to the path, and rmpath to remove the directories from the path. You can have look at the help how to use "path" function in the MATLAB by typing "help path" in the Common Window( see figure below).

fig6.pngfig7.png

Workspace Browser

The MATLAB workspace consists of the set of variables (named arrays) built up during a MATLAB session and stored in memory (see figure below). You add variables to the workspace by using functions, running M-files, and loading saved workspaces.
to view workspace and information about each variable, use the Workspace Browser, or use the function "who" and "whos" , for example, typing

who a*

in the Command Window will show variable names starting with "a".
fig8.png
To delete variables from the workspace, select the variable and select Delete from theEdit menu. Alternatively, use the clear function.
The workspace is not maintained after you end the MATLAB session. To save the workspace to a file that can be read during a later MATLAB session, select Select Workspace As from the File menu. use the save function. This saves the workspace to a binary file called a MAT-file, which has a .mat extension.
There are options for saving to different formats. To read in a MAT-file, select Import Data from the File menu, or use load function.

Array Editor

Double-click a variable in the Workspace browser to see it in the Array Editor (as shown in Figure below). Use the Array Editor to view and edit a visual representation of one- or two-dimensional numeric arrays, strings, and cell arrays of strings that are in the workspace.
fig9.png

Editor/Debugger

Use the Editor/Debugger to create and debug M-files, which are programs you write to run MATLAB functions. The Editor/Debugger provides a graphical user interface for basic text editing, as well as for M-file debugging (see Figure below).

fig10.png

You can use any text editor to create M-files, such as TextPad, and can use preferences (accessible from the desktop File menu) to specify that editor as the default. If you use another editor, you can still use the MATLAB Editor/Debugger for debugging. If you just need to view the contents of an M-file, you can display it in the Command Window by using the type function.


Lab 2: WORKING WITH MATRICES IN MATLAB

Objectives:

Learn about Matrix operations in Matlab
Learn about working with variables, numbers, operators, functions, expressions
Learn about generating matrices, load matrices, create matrices from M-files and concatentation, and delete matrix rows and columns Lear how to use matrices for linear algebra, work with arrays, multivariate data, scalar expansion, and logical subscripting, and use the find function list item

Prerequisites

Students are expected to be familiar with Windows system
Students should have basic knowledge of any text editors

Requirement

Matlab software any version

Remarks

This lab session will be evaluated and it will carry 2 % of your activity credit points.
All the source codes and executables related to this lab session should be saved in a folder named Lab2. If not, zero marks will be given.
Before you leave the lab make sure you show to your instructor that you have saved your file in Lab2 folder.

Matrices in Matlab

In MATLAB, a matrix is a two dimensional array of numbers. The 1-by-1 matrices are called scalars, and matrices with only one row or column are called vectors. MATLAB has other ways of storing both numeric and nonnumeric data, but in the beginning, it is usually best to think of everything as a matrix. The operations in MATLAB are designed to be as natural as possible. Where other programming languages work with numbers one at a time, MATLAB allows you to work with entire matrices quickly and easily.

Entering Matrices

Start MATLAB and follow along with each example. Start by entering a matrix as a list of its elements. You only have to follow a few basic conventions:

  1. Separate the elements of a row with blanks or commas.
  2. Use a semicolon,";", to indicate the end of each row.
  3. Surround the entire list of elements with square brackets [ ].

To enter square matrix, simply type in the Command Window, for example,

A = [1 3 2; 4 2 6; 5 9 8]

And the result should look like (see also figure below):
A =

     1     3     2
     4     2     6
     5     9     8
fig11.png

Once you have entered the matrix, it is automatically remembered in the MATLAB workspace. You can refer to it simply as A. Now that you have A in the workspace, take a look at what makes it so interesting.

sum, transpose, and diag

The first statement to try is

>> sum(A)

MATLAB replies with

ans =

    10    14    16

When you do not specify an output variable, MATLAB uses the variable ans, short for answer, to store the results of a calculation. You have computed a row vector containing the sums of the columns of A.

How about the row sums? MATLAB has a preference for working with the columns of a matrix, so the easiest way to get the row sums is to transpose first the matrix, compute the column sums of the transpose, and then transpose the result. The transpose operation is denoted by an apostrophe or single quote as:

>> A’

It flips a matrix about its main diagonal and it turns a row vector into a column vector:

ans =

     1     4     5
     3     2     9
     2     6     8

Then,

>> sum(A’)’

will produce:

ans =

     6
    12
    22

You can check it easily, that this is the sum along the rows of the matrix A.

The sum of the elements on the main diagonal is obtained with the sum and the diag functions.

>> diag(A)

produces

ans =

     1
     2
     8

and

>> sum(diag(A))

produces

ans =

    11

Subscripts
The element in row i and column j of A is denoted by A(i,j). For example, A(1,2) is the number in the first row and second column. For our square matrix, A(1,2) is 3, which can be verified by typing:

>> A(1,2)
ans =

     3

So to compute the sum of the elements in the first column of A, type:

>> A(1,1) + A(2,1) + A(3,1)

which produces:

ans =

    10

but is not the most elegant way of summing a single column.
It is also possible to refer to the elements of a matrix with a single subscript, A(k). This is the usual way of referencing row and column vectors. But it can also apply to a fully two-dimensional matrix, in which case the array is regarded as one long column vector formed from the columns of the original matrix. So, for our matrix, A(4) is another way of referring to the value
3 stored in A(1,2):

>> A(1,2)

ans =

     3
>> A(4)

ans =

     3

If you try to use the value of an element outside of the matrix, it is an error.

>> A(4,5)

Index exceeds matrix dimensions.

The Colon Operator

The colon, :, is one of the most important MATLAB operators. It occurs in several different forms. The expression
1:5 is a row vector containing the integers from 1 to 5

ans =

     1     2     3     4

5

To obtain non unit spacing, specify an increment. For example,

>> 10:5:30

ans =

    10    15    20    25

30

and

>>0:pi/4:pi

is

0 0.7854 1.5708 2.3562 3.1416

Subscript expressions involving colons refer to portions of a matrix. A(1:k,j) is the first k elements of the jth column of A. So

>>sum(A(1:3,3))

Computes the sum of the third column:

ans =

    16

But there is a better way. The colon by itself refers to all the elements in a row or column of a matrix and the keyword end refers to the last row or column. So

>>sum(A(:,end))

computes the sum of the elements in the last column of A.

ans =

    16

Expressions

Like most other programming languages, MATLAB provides mathematical expressions, but unlike most programming languages, these expressions involve entire matrices. The building blocks of expressions are

  • Variables
  • Numbers
  • Operators
  • Functions

Variables

MATLAB does not require any type declarations or dimension statements. When MATLAB encounters a new variable name, it automatically creates the variable and allocates the appropriate amount of storage. If the variable already exists, MATLAB changes its contents and, if necessary, allocates new storage. For example,

>> i = 25

creates a 1-by-1 matrix named i and stores the value 25 in its single element.

Variable names consist of a letter, followed by any number of letters, digits, or underscores. MATLAB uses only the first 31 characters of a variable name. MATLAB is case sensitive; it distinguishes between uppercase and lowercase letters. A and a are not the same variable. To view the matrix assigned to any variable, simply enter the variable name.

Numbers

MATLAB uses conventional decimal notation, with an optional decimal point and leading plus or minus sign, for numbers. Scientific notation uses the letter e to specify a power-of-ten scale factor. Imaginary numbers use either i or j as a suffix.

Some examples of legal numbers are:

3 -99 0.0001 9.6397238 1.60210e-20 6.02252e23 1i -3.14159j 3e5i

Operators

Expressions use familiar arithmetic operators and precedence rules.

(+) Addition
(-) Subtraction
(*) Multiplication
(/) Division
(\) Left division (described in “Matrices and Linear Algebra” in the MATLAB documentation)
(^) Power
(')Complex conjugate transpose
( ) Specify evaluation order

Functions
MATLAB provides a large number of standard elementary mathematical functions, including abs, sqrt, exp, and sin. Taking the square root or logarithm of a negative number is not an error; the appropriate complex result is produced automatically. MATLAB also provides many more advanced mathematical functions, including Bessel and gamma functions. Most of these functions accept complex arguments. For a list of the elementary mathematical functions, type
help elfun

For a list of more advanced mathematical and matrix functions, type

help specfun
help elmat

Some of the functions, like sqrt and sin, are built in. They are part of the MATLAB core so they are very efficient, but the computational details are not readily accessible. Other functions, like gamma and sinh, are implemented in M-files. You can see the code and even modify it if you want.

Several special functions provide values of useful constants.

pi 3.14159265…
i Imaginary unit, √-1
j Same as i
eps Floating-point relative precision, 2-52
realmin Smallest floating-point number, 2-1022
realmax Largest floating-point number, (2-ε)21023
Inf Infinity
NaN Not-a-number

Infinity is generated by dividing a nonzero value by zero, or by evaluating well defined mathematical expressions that overflow, i.e., exceed realmax. Not-a-number is generated by trying to evaluate expressions like 0/0 or Inf-Inf that do not have well defined mathematical values.

The function names are not reserved. It is possible to overwrite any of them with a new variable, such as

>>eps = 1.e-6

and then use that value in subsequent calculations. The original function can be restored with

>>clear eps

Examples of Expressions

You have already seen several examples of MATLAB expressions. Here are a
few more examples, and the resulting values.

>> r=1+sqrt(2)/5

r =
    1.2828

Working with matrices

Generating Matrices

zeros All zeros
ones All ones
rand Uniformly distributed random elements
randn Normally distributed random elements

>> z=zeros(3,3)

z =

     0     0     0
     0     0     0
     0     0     0
>> T=2*ones(2,3)

T =

     2     2     2
     2     2     2
>> N=rand(1,5)

N =

    0.8147    0.9058    0.1270    0.9134    0.6324
>> G=randn(1,5)

G =

   -1.3077   -0.4336    0.3426    3.5784    2.7694

M-Files

You can create your own matrices using M-files, which are text files containing MATLAB code. Use the MATLAB Editor or another text editor to create a file name it matFile.m.

For example, create a file containing these five lines.

A = [ ...
16.0 3.0 2.0 13.0
5.0 10.0 11.0 8.0
9.0 6.0 7.0 12.0
4.0 15.0 14.0 1.0 ];

Then the statement

matFile

reads the file and creates a variable, A, containing our example matrix.

A =

    16     3     2    13
     5    10    11     8
     9     6     7    12
     4    15    14     1

Linear Algebra

The mathematical operations defined on matrices are the subject of linear algebra. Adding a matrix to its transpose produces a symmetric matrix.

A =

    16     3     2    13
     5    10    11     8
     9     6     7    12
     4    15    14     1
>> A+A'

ans =

    32     8    11    17
     8    20    17    23
    11    17    14    26
    17    23    26     2

The multiplication symbol, *, denotes the matrix multiplication involving inner products between rows and columns. Multiplying the transpose of a matrix by the original matrix also produces a symmetric matrix.

>> A'*A
ans =
   378   212   206   360
   212   370   368   206
   206   368   370   212
   360   206   212   378

The determinant of this particular matrix happens to be zero, indicating that
the matrix is singular.

>> det(A)

ans =

  1.0871e-012

Since the matrix is singular, it does not have an inverse. If you try to compute the inverse with you will get a warning message

>> inv(A)

Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 9.796086e-018.

ans =

  1.0e+015 *

    0.1251    0.3753   -0.3753   -0.1251
   -0.3753   -1.1259    1.1259    0.3753
    0.3753    1.1259   -1.1259   -0.3753
   -0.1251   -0.3753    0.3753    0.1251

The eigenvalues of the magic square are interesting.

>> eig(A)

ans =

   34.0000
    8.0000
    0.0000
   -8.0000

If the square matrix is multiplied by itself with array multiplication A.*A the result is an array containing the squares of the integers from 1 to 16, in an unusual order.

>> A.*A

ans =

   256     9     4   169
    25   100   121    64
    81    36    49   144
    16   225   196     1

Lab 3: WORKING WITH MATRICES IN MATLAB

Objectives:
Learn more about Matrix operations in Matlab
Learn about creating plots in Matlab.

Prerequisites
Students are expected to be familiar with Windows system
Students should have basic knowledge of any text editors

Requirement
Matlab software any version

Remarks
This lab session will be evaluated and it will carry 2 % of your activity credit points.
All the source codes and executables related to this lab session should be saved in a folder named Lab3. If not, zero marks will be given.
Before you leave the lab make sure you show to your instructor that you have saved your file in Lab3 folder.

More about matrices in Matlab

We can reshape matrices in Matlab, for example we will consider reshaping a 3-by-4 matrix A to have dimensions 2-by-6:

>> A = rand(3,4);
>> A

A =

    0.8944    0.9274    0.6183    0.1248
    0.1375    0.9175    0.3433    0.7306
    0.3900    0.7136    0.9360    0.6465
>> B=reshape(A, 2, 6)

B =

    0.8944    0.3900    0.9175    0.6183    0.9360    0.7306
    0.1375    0.9274    0.7136    0.3433    0.1248    0.6465

Transpose A so that the row elements become columns. You can use either the transpose function or the transpose operator (.') to do this:
>> A'

ans =

    0.8944    0.1375    0.3900
    0.9274    0.9175    0.7136
    0.6183    0.3433    0.9360
    0.1248    0.7306    0.6465

Rotate the matrix by 90 degrees:

>> B=rot90(A)

B =

    0.1248    0.7306    0.6465
    0.6183    0.3433    0.9360
    0.9274    0.9175    0.7136
    0.8944    0.1375    0.3900

Flip A in a left-to-right direction:

>> B=fliplr(A)

B =

    0.1248    0.6183    0.9274    0.8944
    0.7306    0.3433    0.9175    0.1375
    0.6465    0.9360    0.7136    0.3900

TF = isfinite(A) returns an array the same size as A containing logical 1 (true) where the elements of the array A are finite and logical 0 (false) where they are infinite or NaN.

>> TF = isfinite(A)

TF =

     1     1     1     1
     1     1     1     1
     1     1     1     1

d = size(X) returns the sizes of each dimension of array A in a vector d with ndims(A) elements. If A is a scalar, which MATLAB software regards as a 1-by-1 array, size(A) returns the vector [1 1]. [m,n] = size(A) returns the size of matrix A in separate variables m and n:

>> [m,n]=size(A)

m =
     3

n =

     4

[V,D] = eig(A) produces matrices of eigenvalues (D) and eigenvectors (V) of a square matrix A, so that A*V = V*D. For example:

>> A=rand(3,3)

A =

    0.8332    0.8352    0.9791
    0.3983    0.3225    0.5493
    0.7498    0.5523    0.3304
>> [V,D]=eig(A)

V =

   -0.7744   -0.4620   -0.3334
   -0.3768    0.8400   -0.4880
   -0.5083   -0.2845    0.8067

D =

    1.8822         0         0
         0   -0.0826         0
         0         0   -0.3135

[row,col,v] = find(X) returns a column or row vector v of the nonzero entries in X, as well as row and column indices. If X is a logical expression, then v is a logical array. Output v contains the non-zero elements of the logical array obtained by evaluating the expression X. For example:

>> A=magic(3)

A =

     8     1     6
     3     5     7
     4     9     2
>> [r,c,v]= find(A>5);
>> r', c', v'

ans =

     1     3     1     2

ans =

     1     2     3     3

ans =

     1     1     1     1

The scalar product of the vectors A and B is computed as:

dot(A,B)

where A and B must be vectors of the same length. When A and B are both column vectors, DOT(A,B) is the same as A'*B:

>> A=rand(5,1)

A =

    0.6947
    0.9727
    0.3278
    0.8378
    0.7391
>> B=rand(5,1)

B =

    0.9542
    0.0319
    0.3569
    0.6627
    0.2815
>> C=dot(A,B)

C =

    1.5741
>> D=A'*B

D =

    1.5741

If A and B are row vectors, then dot(A,B) is equivalent as: A*B’. For example,

>> A=rand(1,5)
A =

    0.2304    0.7111    0.6246    0.5906    0.6604
>> B=rand(1,5)

B =

    0.0476    0.3488    0.4513    0.2409    0.7150
>> C=dot(A,B)

C =

    1.1554
>> D=A*B'

D =

    1.1554

We can also find the inverse of a matrix. You must be careful, however, since the operations are numerical manipulations done on digital computers. In the example, the matrix A is not a full matrix, but matlab's inverse routine will still return a matrix.

>> >> A = [ 1 2 3; 3 4 5; 6 7 8]

A =

     1     2     3
     3     4     5
     6     7     8
>> inv(A)

Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 3.469447e-018. 

ans =

  1.0e+015 *

   -2.7022    4.5036   -1.8014
    5.4043   -9.0072    3.6029
   -2.7022    4.5036   -1.8014

By the way, Matlab is case sensitive. This is another potential source of problems when you start building complicated algorithms.

>> inv(a)
??? Undefined function or variable a.

There are also routines that let you find solutions to equations. For example, if Ax=b and you want to find x, a slow way to find x is to simply invert A and perform a left multiply on both sides (more on that later). It turns out that there are more efficient and more stable methods to do this (L/U decomposition with pivoting, for example). Matlab has special commands that will do this for you.
Before finding the approximations to linear systems, it is important to remember that if A and B are both matrices, then AB is not necessarily equal to BA. To distinguish the difference between solving systems that have a right or left multiply, Matlab uses two different operators, "/" and "\". Examples of their use are given below. It is left as an exercise for you to figure out which one is doing what.

>> v = [1 3 5]'

v =

     1
     3
     5
>> x = A\v

Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 4.565062e-18
x =

   1.0e+15 *

    1.8014
   -3.6029
    1.8014
>> x = B\v

x =

     2
     1
    -1
>> B*x

ans =

     1
     3
     5
>> x1 = v'/B
x1 =

    4.0000   -3.0000    1.0000

>> x1*B

ans =

    1.0000    3.0000    5.0000

Vectors and Matrices and Matlab Functions

If you pass a vector to a predefined math function, it will return a vector of the same size, and each entry is found by performing the specified operation on the corresponding entry of the original vector: (Notice that the second command has a ";" at the end of the line. This tells Matlab that it should not print out the result.)

>> v = [1 2 3]';
>> sin(v)

ans =

    0.8415
    0.9093
    0.1411
>> log(v)

ans =

         0
    0.6931
    1.0986

The ability to work with these vector functions is one of the advantages of Matlab. Now complex operations can be defined that can be done quickly and easily. In the following example a very large vector is defined and can be easily manipulated.

>> x = [0:0.1:100]

x =

  Columns 1 through 7 

         0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000

    [stuff deleted]

  Columns 995 through 1001 

   99.4000   99.5000   99.6000   99.7000   99.8000   99.9000  100.0000
>> y = sin(x).*x./(1+cos(x));

Creating a Plot
The plot function has different forms, depending on the input arguments. If y is a vector, plot(y) produces a piecewise linear graph of the elements of y versus the index of the elements of y. If you specify two vectors as arguments, plot(x,y) produces a graph of y versus x.

For example, these statements use the colon operator to create a vector of x values ranging from zero to compute the cosine of these values, and plot the result.

>> x=0:0.1:(2*pi);
>> y=cos(x);
>> plot(x,y)

We can also label the axes and add a title as the following:

>> xlabel('values of x');
>> ylabel('cos(x)');
>> title('cosine function plot')
The plot should look like one in figure below.
figure12.png
You can also specify the font size as follows:
>> xlabel('values of x','FontSize',12);
>> ylabel('cos(x)', 'FontSize',12);
>> title('cosine function plot','FontSize',12)

Multiple Data Sets in One Graph
Multiple x-y pair arguments create multiple graphs with a single call to plot. MATLAB automatically cycles through a predefined (but user settable) list of colors to allow discrimination among sets of data. For example, these statements plot three related functions of x, each curve in a separate distinguishing color.

>> y2 = cos(5*x);
>> y3 = cos(2.5*x);
>> plot(x,y,x,y2,x,y3)

The legend command provides an easy way to identify the individual plots (see also Figure 2).

>> legend('cos(x)','cos(5*x)','cos(2.5*x)')

figure13.png

Specifying Line Styles and Colors
It is possible to specify color, line styles, and markers (such as plus signs or circles) when you plot your data using the plot command.

plot(x,y,'color_style_marker')

color_style_marker is a string containing from one to four characters (enclosed in single quotation marks) constructed from a color, a line style, and a marker type:

Color strings are 'c', 'm', 'y', 'r', 'g', 'b', 'w', and 'k'. These correspond to cyan, magenta, yellow, red, green, blue, white, and black.

Linestyle strings are '-' for solid, '—' for dashed, ':' for dotted, '-.' For dash-dot. Omit the linestyle for no line.
The marker types are '+', 'o', '*', and 'x' and the filled marker types are's' for square, 'd' for diamond, '^' for up triangle, 'v' for down triangle, '>' for right triangle, '<' for left triangle, 'p' for pentagram, 'h' for hexagram, and none for no marker.

Plotting Lines and Markers
If you specify a marker type but not a linestyle, MATLAB draws only the marker. For example,

plot(x,y,'ko')

plots black circle at each data point, but does not connect the markers with a line. The statement

plot(x,y,'b-o')

plots a blue solid line and places circle sign markers at each data point. You may want to use fewer data points to plot the markers than you use to plot the lines. This example plots the data twice using a different number of points for the solid line and marker plots (see Figure below).

x1 = 0:pi/100:2*pi;
y1=cos(x1);
x2 = 0:pi/10:2*pi;
y2=cos(x2);
plot(x1,y1,'r-',x2,y2,'ro')

figure14.png

Adding Plots to an Existing Graph
The hold command enables you to add more than one plot to an existing graph. When you type

hold on

after each plot command MATLAB does not replace the existing graph, but it adds the new data to the current graph, rescaling the axes if needed.

For example the following Matlab commands will create the same plot as in previous figure.

x1 = 0:pi/100:2*pi;
y1=cos(x1);
plot(x1, y1, ‘r-‘);
hold on;
x2 = 0:pi/10:2*pi;
y2=cos(x2);
plot(x2,y2,'ro')
hold on;

Lab4: WORKING WITH PLOTS AND CODING IN MATLAB

Objectives:
Learn more about creating plots in Matlab.
Learn about subplot.
Learn about coding in Matlab.

Prerequisites
Students are expected to be familiar with Windows system
Students should have basic knowledge of any text editors
Requirement

Matlab software any version
Remarks
This lab session will be evaluated and it will carry 2 % of your activity credit points.
All the source codes and executables related to this lab session should be saved in a folder named Lab4. If not, zero marks will be given.
Before you leave the lab make sure you show to your instructor that you have saved your file in Lab4 folder.

More about plots in Matlab

Figure Windows
If there are no figure windows already on the screen, plotting functions automatically open a new figure window. If a figure window exists, MATLAB uses that window for graphics output. If there are multiple figure windows open, MATLAB targets the one that is designated the “current figure” (the last figure used or clicked in).
To make an existing figure window the current figure, you can click the mouse while the pointer is in that window or you can type

figure(n)

where n is the number in the figure title bar. The results of subsequent graphics commands are displayed in this window.

To open a new figure window and make it the current figure, type

figure

Multiple Plots in One Figure
The subplot command enables one to display multiple plots in the same figure or print them on the same piece of paper. Typing

subplot(m,n,p)

partitions the figure window into an m-by-n matrix of small subplots and selects the pth subplot for the current plot. The axes are counted along the top row of the figure window, then the second row, etc. For example,

income = [3.2 4.1 5.0 5.6];
outgo = [2.5 4.0 3.35 4.9];
subplot(2,1,1); plot(income)
title('Income')
subplot(2,1,2); plot(outgo)
title('Outgo')
plots income on the top half of the window and outgo on the bottom half as shown in figure below.
fig15.png

Another example is given below creating 4 subplots (see also Figure 2):

t= 0:pi/10:2*pi;
[X,Y,Z] = cylinder(4*cos(t));
subplot(2,2,1); mesh(X)
subplot(2,2,2); mesh(Y)
subplot(2,2,3); mesh(Z)
subplot(2,2,4); mesh(X,Y,Z)
fig16.png

Controlling the Axes
The axis command supports a number of options for setting the scaling, orientation, and aspect ratio of plots.

Setting Axis Limits
The axis command enables you to specify the axes limits, for example:

axis([xmin xmax ymin ymax])

sets the limits for the x- and y-axis of the current axes. For three-dimensional plot, we may use:

axis([xmin xmax ymin ymax zmin zmax])

While the command

axis auto

enables MATLAB automatic limit selection.

Setting Axis Aspect Ratio

axis also enables you to specify a number of predefined modes. For example,

axis square

makes the x-axes and y-axes the same length.

axis equal

makes the individual tick mark increments on the x- and y-axes the same length. This means

plot(exp(i*[0:pi/10:2*pi]))

followed by either axis square or axis equal turns the oval into a proper circle.

axis auto normal

returns the axis scaling to its default, automatic mode.

Setting Axis Visibility

You can use the axis command to make the axis visible or invisible.

axis on

makes the axis visible. This is the default.

axis off

makes the axis invisible.

Setting Grid Lines

The grid command toggles grid lines on and off. The statement

grid on

turns the grid lines on and

grid off

turns them back off again.

Axis Labels and Titles
The xlabel, ylabel, and zlabel commands add x-, y-, and z-axis labels, respectively. The title command adds a title at the top of the figure and the text function inserts text anywhere in the figure. A subset of TeX notation produces Greek letters. You can also set these options interactively.

For example,

t = -pi:pi/100:pi;
y = sin(t).*cos(t);
plot(t,y)
axis([-pi pi -1 1])
xlabel('-\pi \leq {\itt} \leq \pi')
ylabel('sin(t)*cos(t)')
title('Graph of the sine \times cosine functions')
text(-2.5,0.6,'{\itMax}')
text(-0.9,-0.6,'{\itMin}')
text(0.9,0.6,'{\itMax}')
text(2.5,-0.6,'{\itMin}')
fig17.png

You can save the figure in different formats. To do this, select Save from the File menu. To save it using a graphics format, such as JPEG, PNG, etc., for use with other applications, select Export from the File menu, then graph shown in figure will appear. Select Export to save the file using the desired format.

fig18.png

Programming with Matlab

Flow Control

MATLAB has several flow control constructs:
• “if”
• “switch and case”
• “for”
• while
• “continue”
• “break”

if

The if statement evaluates a logical expression and executes a group of statements when the expression is true. The optional elseif and else keywords provide for the execution of alternate groups of statements. An end keyword, which matches the if, terminates the last group of statements. The groups of statements are delineated by the four keywords—no braces or brackets are involved.

The MATLAB algorithm for generating a random square n-by-n involves three different cases: when n is odd, when n is even but not divisible by 4, or when n is divisible by 4. This is described by

if rem(n,2) ~= 0
M = randn(n,n);
elseif rem(n,4) ~= 0
M = ones(n,n)
else
M = zeros(n,n)
End

In this example, the three cases are mutually exclusive, but if they weren’t, the first true condition would be executed. REM function is the remainder after division.

switch and case

The switch statement executes groups of statements based on the value of a variable or expression. The keywords case and otherwise delineate the groups. Only the first matching case is executed. There must always be an end to match the switch.

The previous example can also be written as:

switch (rem(n,4)==0) + (rem(n,2)==0)
case 0
M = randn(n,n)
case 1
M = ones(n,n)
case 2
M = zeros(n,n)
otherwise
error('This is impossible')
end

for

The for loop repeats a group of statements a fixed, predetermined number of times. A matching end delineates the statements.

for n = 3:32
    r(n) = randn(1,1);
end
r

The semicolon terminating the inner statement suppresses repeated printing, and the r after the loop displays the final result. It is a good idea to indent the loops for readability, especially when they are nested. For this, use Text →Smart Indent.

for i = 1:m
     for j = 1:n 
          H(i,j) = 1/(i+j);
     end
end

while

The while loop repeats a group of statements an indefinite number of times under control of a logical condition. A matching end delineates the statements.

i=0;
n=20;
s=0;
while i < n
    i=i+1;
    s = s+i^2;
end
s

continue

The continue statement passes control to the next iteration of the for or while loop in which it appears, skipping any remaining statements in the body of the loop. In nested loops, continue passes control to the next iteration of the for or while loop enclosing it.

The example below shows a continue loop that counts the lines of code in the file, magic.m, skipping all blank lines and comments. A continue statement is used to advance to the next line in magic.m without incrementing the count whenever a blank line or comment line is encountered.

fid = fopen('splots6.m','r');
count = 0;
while ~feof(fid)
    line = fgetl(fid);
    if isempty(line) | strncmp(line,'%',1)
        continue
    end
    count = count + 1;
end
disp(sprintf('%d lines',count));

where ``splots6.m” looks as follows:

i=0;
n=20;
s=0;
while i < n
    i=i+1;
    s = s+i^2;
end
s

break

The break statement lets you exit early from a for or while loop. In nested loops, break exits from the innermost loop only. Here is an improvement on the example from the previous section. Why is this use of break a good idea?

i=0;
n=20;
s=0;
while i < n
    i=i+1;
    s = s+i^2;
    if s>= 100 
        break;
    end
end
s

Lab 5: SYSTEM MODELING AND NUMERICAL INTEGRATION IN MATLAB

Objectives:

Learn more about modelling a simple system in matlab.
Learn about Euler method of numerical integration in matlab.

Prerequisites

Students are expected to be familiar with Windows system
Students should have basic knowledge of any text editors

Requirement

Matlab software any version

Remarks

This lab session will be evaluated and it will carry 2 % of your activity credit points.
All the source codes and executables related to this lab session should be saved in a folder named Lab5. If not, zero marks will be given.
Before you leave the lab make sure you show to your instructor that you have saved your file in Lab5 folder.

Example 1: Conveyor System

Consider a factory conveyor system in which boxes arrive at the rate of one box each Δt. Each box has a weight of w1, w2 , w3 . However, there are twice as many w1 boxes and w3 boxes as w2 boxes. How do we model and simulate this?

fig19.png

The Probability mass distribution is:

(1)
\begin{equation} p_1 = P(W=w_1)=0.4, p_2 = P(W=w_2)=0.2, p_3 = P(W=w_3)=0.4. \end{equation}

where W is a weight random variable that can take on one of the three discrete values. The notation P(W=w) gives the probability that the random variable W is w.
The set {w1 w2 w3 } is called the sample space of , and is the set of all possible weights. According to the description, these boxes arrive every Δt seconds, so t = Δt k gives the continuous time measured in successive k-values, assuming the initial time is zero. However, how do we describe the system output?

An interesting problem would be the weight W(t) of the boxes as they arrive. W(t) is a non-deterministic variable, and we can only hope to simulate the behavior of this variable. If we denote with W(k) the weight of the k-th event. This can be accomplished by using the random (RAND) function, which is a random number generator that provides uniformly random distributed variates such that

(2)
\begin{equation} 0≤rand≤1 \end{equation}
function lab5_ex1 
clear all; 
N=100; 
w=zeros(1,N); 
w1=1; w2=2; w3=3; 
Nw=zeros(1,3); 
P=zeros(1,3); 
for k=1:N
     r=10*rand(1,1);
     if r < 4 
            w(k) = w1; 
            Nw(1)=Nw(1)+1; 
     elseif r >= 4 && r < 6
              w(k) = w2; 
              Nw(2)=Nw(2)+1; 
     else 
             w(k) = w3;
             Nw(3)=Nw(3)+1;
     end 
end 
P=Nw/N; 
Pw=[0.4 0.2 0.4];
figure(1); 
hold on;
subplot(2,1,1); bar(P,'grouped','r'); 
subplot(2,1,2); bar(Pw,'grouped','b'); 
colormap jet; 

figure(2); 
hold on; 
stairs(cumsum(w), 'k-');

Euler Numerical integration method:
A matlab code for numerical integration using Euler 1st order method is shown below:

function [wi, ti] = euler ( RHS, t0, x0, tf, N ) 

%EULER approximate the solution of the initial value problem 
% 
%         x'(t) = RHS( t, x ), x(t0) = x0 
% 
%    using Euler's method - this routine will work for a system
%    of first-order equations as well as for a single equation 
% 
% calling sequences: 
%        [wi, ti] = euler ( RHS, t0, x0, tf, N ) 
%        euler ( RHS, t0, x0, tf, N ) 
% 
% inputs: 
%       RHS   string containing name of m-file defining the 
%           right-hand side of the differential equation; the 
%           m-file must take two inputs - first, the value of 
%           the independent variable; second, the value of the 
%           dependent variable
%      t0   initial value of the independent variable 
%      x0   initial value of the dependent variable(s) 
%          if solving a system of equations, this should be a 
%          row vector containing all initial values 
%     tf   final value of the independent variable 
%     N    number of uniformly sized time steps to be taken to 
%          advance the solution from t = t0 to t = tf 
% 
% output:
%     wi    vector / matrix containing values of the approximate 
%          solution to the differential equation 
%     ti     vector containing the values of the independent 
%          variable at which an approximate solution has been 
%          obtained 
%
neqn = length ( x0 ); 
ti = linspace ( t0, tf, N+1 ); 
wi = [ zeros( neqn, N+1 ) ]; 
wi(1:neqn, 1) = x0'; 

h = ( tf - t0 ) / N; 
for i = 1:N 
   x0 = x0 + h * feval ( RHS, t0, x0 ); 
           t0 = t0 + h;

   wi(1:neqn,i+1) = x0'; 
end;

Example

Consider the following differential equation:

(3)
\begin{align} \frac{dx}{dt} = x^2*t \end{align}

With initial conditions:

(4)
\begin{equation} t_0=1;x(t_0)=3 \end{equation}

A matlab function for returning a value of $\frac{dx}{dt}$ is given in the following:

function f = example(t,x) 
f=t*(x^2);
end

A matlab testing Euler method for different values of the time step is given below:

function Eulerdemo 

x0=3; 
t0=1; 
tf=1.3; 

h=0.02; 
N=floor((tf-t0)/h) 
 [x4, t4]=euler('example', t0, x0, tf, N); 

xt = 6./(5-3*t4.^2); 

plot(t4,x4, 'r-x', 'LineWidth',2,...
                   'MarkerEdgeColor','r',... 
                   'MarkerFaceColor','r',... 
                   'MarkerSize',10); 
hold on; 

plot(t4, xt, 'k-', 'LineWidth',2); 
hold on; 

h=0.04; 
N=floor((tf-t0)/h) 
 [x4, t4]=euler('example', t0, x0, tf, N); 

plot(t4,x4, 'b-o', 'LineWidth',2,... 
                    'MarkerEdgeColor','r',... 
                    'MarkerFaceColor','r',... 
                    'MarkerSize',10);
hold on; 

h=0.08; 
N=floor((tf-t0)/h) 
 [x4, t4]=euler('example', t0, x0, tf, N); 

plot(t4,x4, 'g-s', 'LineWidth',2,...
                   'MarkerEdgeColor','r',... 
                   'MarkerFaceColor','r',... 
                   'MarkerSize',10); 
hold on; 

xlim([1 1.25]);

Graphically results are presented in figure:

fig20.png

The analytical solution of the above differential equation is:

(5)
\begin{align} x(t)=\frac{6}{5-3t^2} \end{align}

+ Lab 6: TAYLOR AND RUNGE-KUTTA NUMERICAL INTEGRATORS IN MATLAB


Objectives:

  • Learn more about numerical integration in matlab.
  • Comparison of numerical integrators in matlab.

Prerequisites

  • Students are expected to be familiar with Windows system
  • Students should have basic knowledge of any text editors

Requirement

Matlab software any version

Remarks

This lab session will be evaluated and it will carry 2 % of your activity credit points.
All the source codes and executables related to this lab session should be saved in a
folder named Lab6. If not, zero marks will be given.
Before you leave the lab make sure you show to your instructor that you have saved
your file in Lab6 folder.

Example 1:

Consider again the following differential equation:

(6)
\begin{align} \frac{dx}{dt}=x^2t \end{align}

with initial condition:

(7)
\begin{equation} x(1)=3. \end{equation}

A matlab function for returning a value of $\frac{dx}{dt}$ is given in the following:

function f = example(t,x) 
 f=t*(x^2); 
end

Euler Numerical integration method:

A matlab code for numerical integration using Euler 1st order method is shown below:

function [wi, ti] = euler ( RHS, t0, x0, tf, N ) 

%EULER approximate the solution of the initial value problem 
% 
% x'(t) = RHS( t, x ), x(t0) = x0 
% 
% using Euler's method - this routine will work for a system 
% of first-order equations as well as for a single equation 
% 
% calling sequences: 
% [wi, ti] = euler ( RHS, t0, x0, tf, N ) 
% euler ( RHS, t0, x0, tf, N ) 
% 
% inputs: 
% RHS string containing name of m-file defining the 
% right-hand side of the differential equation; the 
% m-file must take two inputs - first, the value of 
% the independent variable; second, the value of the 
% dependent variable 
% t0 initial value of the independent variable 
% x0 initial value of the dependent variable(s) 
% if solving a system of equations, this should be a 
% row vector containing all initial values 
% tf final value of the independent variable 
% N number of uniformly sized time steps to be taken to 
% advance the solution from t = t0 to t = tf 
% 
% output: 
% wi vector / matrix containing values of the approximate 
% solution to the differential equation 
% ti vector containing the values of the independent 
% variable at which an approximate solution has been 
% obtained 
% 

neqn = length ( x0 ); 
ti = linspace ( t0, tf, N+1 ); 
wi = [ zeros( neqn, N+1 ) ]; 
wi(1:neqn, 1) = x0'; 

h = ( tf - t0 ) / N; 

for i = 1:N 
 x0 = x0 + h * feval ( RHS, t0, x0 ); 
 t0 = t0 + h; 

 wi(1:neqn,i+1) = x0'; 
end;
Taylor 2nd order numerical method in matlab: 

function [wi, ti] = taylor2nd ( RHS, t0, x0, tf, N ) 
%TAYLOR2ND approximate the solution of the initial value problem 
% S
% x'(t) = RHS( t, x ), x(t0) = x0 
% 
% using the second-order Taylor method - this routine will 
% work for a system of first-order equations as well as for 
% a single equation 
% 
% calling sequences: 
% [wi, ti] = taylor2nd (RHS, t0, x0, tf, N ) 
% taylor2nd ( RHS, t0, x0, tf, N ) 
% 
% inputs: 
% RHS string containing name of m-file defining the 
% right-hand side of the differential equation and 
% its derivative with respect to the independent 
% variable; 
% the prototype for the m-file should be 
% [f, ft] = RHS ( t, x ) 
% where f is the value of the right-hand side function 
% and ft is the value of the derivative with respect 
% to the independent variable 
% t0 initial value of the independent variable 
% x0 initial value of the dependent variable(s) 
% if solving a system of equations, this should be a 
% row vector containing all initial values 
% tf final value of the independent variable 
% N number of uniformly sized time steps to be taken to 
% advance the solution from t = t0 to t = tf 
% 
% output: 
% wi vector / matrix containing values of the approximate
solution to the differential equation 
% ti vector containing the values of the independent 
% variable at which an approximate solution has been 
% obtained 
% 

neqn = length ( x0 ); 
ti = linspace ( t0, tf, N+1 ); 
wi = [ zeros( neqn, N+1 ) ]; 
wi(1:neqn, 1) = x0'; 

h = ( tf - t0 ) / N; 

for i = 1:N 
 [f ft] = feval ( RHS, t0, x0 ); 
 x0 = x0 + h * f + h^2 * ft / 2; 
 t0 = t0 + h; 

 wi(1:neqn,i+1) = x0'; 
end;

A matlab testing Euler method for different values of the time step is given below:

function Eulerdemo 
clear all;
x0=3; 
t0=1; 
tf=1.3;  
h=0.02; 
N=floor((tf-t0)/h) 
[x4, t4]=euler('example', t0, x0, tf, N); 
xt = 6./(5-3*t4.^2); 
plot(t4,x4, 'r-x', 'LineWidth',2,... 
                      'MarkerEdgeColor','r',... 
                      'MarkerFaceColor','r',... 
                      'MarkerSize',10); 
hold on; 
plot(t4, xt, 'k-', 'LineWidth',2); 
hold on;
h=0.04; 
N=floor((tf-t0)/h) 
[x4, t4]=euler('example', t0, x0, tf, N); 
plot(t4,x4, 'b-o', 'LineWidth',2,... 
                       'MarkerEdgeColor','r',... 
                       'MarkerFaceColor','r',... 
                       'MarkerSize',10); 
hold on; 
h=0.08; 
N=floor((tf-t0)/h) 
[x4, t4]=euler('example', t0, x0, tf, N); 
hold on; 

xlim([1 1.25]);

Runge-Kutta 4th order:

function [wi, ti] = rk4 ( RHS, t0, x0, tf, N ) 
%RK4 approximate the solution of the initial value problem 
% 
% x'(t) = RHS( t, x ), x(t0) = x0 
% 
% using the classical fourth-order Runge-Kutta method - this 
% routine will work for a system of first-order equations as 
% well as for a single equation 
% 
% calling sequences:

plot(t4,x4, 'g-s', 'LineWidth',2,... 
                       'MarkerEdgeColor','r',... 
                       'MarkerFaceColor','r',... 
                       'MarkerSize',10);
% [wi, ti] = rk4 ( RHS, t0, x0, tf, N ) 
% rk4 ( RHS, t0, x0, tf, N ) 
% 
% inputs: 
% RHS string containing name of m-file defining the 
% right-hand side of the differential equation; 
the 
% m-file must take two inputs - first, the value 
of 
% the independent variable; second, the value of 
the 
% dependent variable 
% t0 initial value of the independent variable 
% x0 initial value of the dependent variable(s) 
% if solving a system of equations, this should 
be a 
% row vector containing all initial values 
% tf final value of the independent variable 
% N number of uniformly sized time steps to be 
taken to 
% advance the solution from t = t0 to t = tf 
% 
% output: 
% wi vector / matrix containing values of the 
approximate 
% solution to the differential equation 
% ti vector containing the values of the independent 
% variable at which an approximate solution has 
been 
% obtained 
% 

neqn = length ( x0 ); 
ti = linspace ( t0, tf, N+1 ); 
wi = [ zeros( neqn, N+1 ) ]; 
wi(1:neqn, 1) = x0'; 

h = ( tf - t0 ) / N; 

for i = 1:N 
 k1 = h * feval ( RHS, t0, x0 ); 
 k2 = h * feval ( RHS, t0 + h/2, x0 + k1/2 ); 
 k3 = h * feval ( RHS, t0 + h/2, x0 + k2/2 ); 
 k4 = h * feval ( RHS, t0 + h, x0 + k3 ); 
 x0 = x0 + ( k1 + 2*k2 + 2*k3 + k4 ) / 6; 
 t0 = t0 + h; 

 wi(1:neqn,i+1) = x0'; 
end;

Adaptive Runge-Kutta method:

function [wi, ti, count] = rkf45 ( RHS, t0, x0, tf, parms ) 

%RKF45 approximate the solution of the initial value problem 
% 
% x'(t) = RHS( t, x ), x(t0) = x0 
% 
% using Runge-Kutta-Fehlberg 4th order - 5th order method - 
% this routine will work for a system of first-order 
    equations 
% as well as for a single equation 
% 
% calling sequences: 
% [wi, ti] = rkf45 ( RHS, t0, x0, tf, parms ) 
% rkf45 ( RHS, t0, x0, tf, parms ) 
% 
% inputs: 
% RHS string containing name of m-file defining the 
% right-hand side of the differential equation; 
   the 
% m-file must take two inputs - first, the value 
    of 
% the independent variable; second, the value of 
    the 
% dependent variable 
% t0 initial value of the independent variable 
% x0 initial value of the dependent variable(s) 
% if solving a system of equations, this should 
    be a 
% row vector containing all initial values 
% tf final value of the independent variable 
% parms three component vector of paramter values 
% parm(1) minimum allowed step size 
% parm(2) maximum allowed step size 
% parm(3) absolute error tolerance 
% 
% output: 
% wi vector / matrix containing values of the 
   approximate 
% solution to the differential equation 
% ti vector containing the values of the independent 
% variable at which an approximate solution has 
    been 
% obtained 
% count number of function evaluations used in 
   advancing the 
% solution from t = t0 to t = tf 
% 

neqn = length ( x0 ); 
hmin = parms(1); 
hmax = parms(2); 
TOL = parms(3); 

ti(1) = t0; 
wi(1:neqn, 1) = x0'; 
count = 0; 
h = hmax; 
i = 2; 

while ( t0 < tf ) 
 k1 = h * feval ( RHS, t0, x0 ); 
 k2 = h * feval ( RHS, t0 + h/4, x0 + k1/4 ); 
 k3 = h * feval ( RHS, t0 + 3*h/8, x0 + 3*k1/32 + 9*k2/32 ); 
 k4 = h * feval ( RHS, t0 + 12*h/13, x0 + 1932*k1/2197 -
7200*k2/2197 + 7296*k3/2197 ); 
 k5 = h * feval ( RHS, t0 + h, x0 + 439*k1/216 - 8*k2 + 
3680*k3/513 - 845*k4/4104 ); 
 k6 = h * feval ( RHS, t0 + h/2, x0 - 8*k1/27 + 2*k2 - 
3544*k3/2565 + 1859*k4/4104 - 11*k5/40 ); 

 R = max ( abs ( k1/360 - 128*k3/4275 - 2197*k4/75240 + k5/50 + 
2*k6/55 ) / h ); 
 q = 0.84 * ( TOL / R ) ^ (1/4); 
 count = count + 6; 

 if ( R < TOL ) 
 x0 = x0 + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50 
+ 2*k6/55; 
% x0 = x0 + 25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5; 
 t0 = t0 + h; 

 ti(i) = t0; 
 wi(1:neqn, i) = x0'; 
 i = i + 1; 
 end; 

 h = min ( max ( q, 0.1 ), 4.0 ) * h; 
 if ( h > hmax ) h = hmax; end; 
 if ( t0 + h > tf ) 
 h = tf - t0; 
 elseif ( h < hmin ) 
 disp ( 'Solution requires step size smaller than minimum' ); 
 return; 
 end; 
end;

A matlab code for testing all the above mentioned methods:

function testRKM 

x0=3; 
t0=1; 
tf=1.3; 
h=0.05; 
N=floor((tf-t0)/h); 
[x4, t4]=euler('example22', t0, x0, tf, N); 
xt = 6./(5-3*(0:0.01:1.25).^2); 
plot((0:0.01:1.25), xt, 'k-', 'LineWidth',2); 
hold on; 
plot(t4,x4, 'r-x', 'LineWidth',2,... 
                      'MarkerEdgeColor','r',... 
                      'MarkerFaceColor','r',... 
                      'MarkerSize',10); 
hold on; 
h=0.05; 
N=floor((tf-t0)/h); 
[x4, t4]=taylor2nd('exampletaylor22', t0, x0, tf, N); 
plot(t4,x4, 'b-o', 'LineWidth',2,... 
                       'MarkerEdgeColor','b',... 
                       'MarkerFaceColor','b',... 
                       'MarkerSize',10); 
hold on; 
h=0.05; 
N=floor((tf-t0)/h);
[x4, t4]=rk4('example22', t0, x0, tf, N); 
plot(t4,x4, 'g-s', 'LineWidth',2,... 
                       'MarkerEdgeColor','g',... 
                       'MarkerFaceColor','g',... 
                       'MarkerSize',10); 
hold on; 
hmin=0.000001; 
hmax = 0.5; 
tol=1.0e-08; 
[x45, t45, c]=rkf45('example22', t0, x0, tf, [hmin hmax tol]); 
plot(t45,x45, 'm o', 'LineWidth',2,... 
                          'MarkerEdgeColor','m',... 
                          'MarkerFaceColor','m',... 
                          'MarkerSize',10); 
hold on; 

xlim([1 1.251]);

Graphically results are presented in figure:

fig21.png

Figure: Comparison of different numerical integration methods for the time step of 0.05.

Lab 7:

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License