CE 206: Engineering Computation S i C t ti Sessional l

1.50 Credits, 3hrs/week

Dr. Tanvir Ahmed Assistant Professor

Course Outline
Introduction to hi-level computational programming tools
- MATLAB, Mathematica etc.

Application to numerical analysis

- B i matrix computations Basic ti t ti - Solving system of linear equations - Solving non-linear equations - Solving differential Equations - Interpolation and curve-fitting - Numerical differentiation - Numerical integration

Application to engineering problems

- Solving mechanics problems - Numerical solution of equation of motion etc. q
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

References and Resources

The MATLAB help file documentation
- Provided with the MATLAB software - www.mathworks.com

Course materials in my website y

- Lecture slides for CE206 - An introduction to MATLAB David Griffith - Getting Started with MATLAB Mathworks Inc Getting MATLAB Inc.

You will also need:

- Your CE205 Course references/Lecture slides

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Course Layout
Problem sets/Assignments
- Approximately one per week - Individual submission in printed form - Include codes and figures

Class grading

- Approximately one per week - Solving a given problem in a given amount of time - Print out the codes and figures and keep it in file

Final Project/Assignment
- Individual or group submission

Requirements
- Basic familiarity with programming (e.g. CE204) - Knowledge on linear algebra, differential equations d i l h d and numerical methods (CE 205)
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Introduction to MATLAB
MATLAB (MATrix LABoratory) is a fully-functional programming language Original code written by Cleve Moler of UNM in the 1970s, later released as a commercial package b M th l d i l k by Mathworks, I k Inc. g y Originally intended for interactive numerical computations I t f ith itt i Interfaces with programs written in other languages(C++, Fortran) Widely used in academia, research institutions and industry worldwide

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

MATLAB basics
MATLAB can be thought of as a super-powerful graphing calculator with many more buttons (i.e. built-in functions) You can build up your own set of functions suited for a particular operation It is an interpreted programming language
-commands are executed line by line

It is capable of graphically representing computational results simultaneously.

- Not possible in C++, Fortran

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

MATLAB as a calculator
Basic arithmetic operators (+, , *, /) used in conjunction with brackets ( )
-5/(4.8+5.32)^2 ans = -0.0488 (3+4i)*(3-4i) ans = 25 cos(pi/2) ans = 6.1230e-017 exp(acos(0.3)) ans = 3 5470 3.5470

5 (4.8 + 5.32) 2

Built in Built-in functions sin() asin() i () i () cos() acos() tan() atan() log() log10() exp() sqrt() abs() round()

(3 + 4i )(3 4i )

cos( / 2) (

cos 1 ( 0.3)

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Numbers and Formats

All computations in MATLAB are done in double precision (15 significant figures)
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Variables
Variable names and their types do not have to be declared in MATLAB. If a statement i t t t t is terminated with a semicolon ( ; ) th results i t d ith i l ), the lt are suppressed letters, Variable names must start with a letter followed by letters digits, and underscores. p The name of variable is not accepted if it is reserved word.

Example
>> x=-13; y = 5*x z = x^2+y x= 13; 5 x, x 2+y y = -65 z = 104

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Variables
Commands involving variables
who: lists the names of the defined variables whos: lists the names and sizes of defined variables clear: clears all variables clear name: clears the variable name clc: clears the command window

Avoid using
ans: default variable name for the result. pi: = 3.1415926 eps: = 2.2204e-016, smallest value by which two numbers 2 2204e 016 can differ inf: , infinity NAN or nan: not-a-number
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Vectors
- Entries within a row are separated by spaces or commas. -Rows are separated by semicolons. - Vector properties using size( ) and length( )
>> a = [1 2 3 4 5 ] a = 1 2 3 >> b = [2;4;6;8;10] b = 2 4 6 8 10 >> si e(a) size(a) ans = 1 >> length(a) ans = 5

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Other methods of vector creation

The colon operator: generally a : b : c produces a vector of entries starting with the value a, incrementing by the value b until it gets t c til t to
>> 3:7 ans = 3 4 5 6 7

>> 0.32:0.1:0.6 ans = 0.3200 0.4200 0.5200 >> -1.4:-0.3:-2 ans = -1.4000 -1.7000 -2.0000

linspace (a,b,n) generates n + 1 equispaced points between a and b inclusive. b, inclusive

CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Graphics: plotting functions

y = sin 3x
for 0 x 1 f
>> x = linspace (0,1,11); >> y = sin(3*pi*x); >> plot(x y) plot(x,y) Increasing the number of elements in x

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Multiplots, titles, labels, linestyles and colors

>> >> >> >> >> plot(x,y,'w-',x,cos(3*pi*x),'g--') legend('Sin curve','Cos curve') title('Multi-plot ') title( Multi plot ) xlabel('x axis'), ylabel('y axis') grid

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Subplot example

>> >> >> >> >> >> >> >>

subplot(221), plot(x,y) b l t(221) l t( ) xlabel('x'),ylabel('sin 3 pi x') subplot(222), plot(x,cos(3*pi*x)) xlabel('x'),ylabel('cos xlabel('x') ylabel('cos 3 pi x') subplot(223), plot(x,sin(6*pi*x)) xlabel('x'),ylabel('sin 6 pi x') subplot(224) plot(x,cos(6*pi*x)) subplot(224), plot(x cos(6*pi*x)) xlabel('x'),ylabel('cos 6 pi x')

subplot(m, n, p) subplot splits the figure window into an mxn array of small axes and makes the pth one active Note - the first subplot active. is at the top left, then the numbering continues across the row
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Subplot example

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Formatted text on plots

plot th first 100 terms in th sequence { n} given by xn = l t the fi t t i the {x i b [1 + 1/n]n and then graph the function (x) = x3 sin2(3x) on the interval -1 x 1
>> set(0,'Defaultaxesfontsize',16); >> n = 1:100; x = (1+1 /n) ^n; (1+1./n). n; >> subplot (211) >> plot(n,x,'.',[0 max(n)],exp(1)*[1 1],... '--','markersize',8) ' ' ' k i ' 8) >> title('x_n = (1+1/n)^n','fontsize',12) >> xlabel('n'), ylabel('x_n') >> legend('x_n','y = e^1 = 2.71828...',4) >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> subplot (212) >> x = -2:.02:2; y = x.^3.*sin(3*pi*x).^2; >> plot(x,y,'linewidth',2) >> legend('y = x^3sin^2 3\pi x',4) g ( y \p , ) >> xlabel('x)
CE 206: Engg. Computation Sessional

Default font changed

Markersize changed Subscript/ superscript LineWidth changed Latin characters

Dr. Tanvir Ahmed

Formatted text on plots

plot th first 100 terms in th sequence { n} given by xn = l t the fi t t i the {x i b [1 + 1/n]n and then graph the function (x) = x3 sin2(3x) on the interval -1 x 1

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Adding additional plots

hold h ld on and h ld off d hold ff hold on tells MATLAB to keep the current data plotted and add y plot graph. the results of any further p commands to the g p Hold off releases the hold on the figure

x = 0:.1:2*pi; 0:.1:2 pi; y = sin(x); p y plot(x,y,'b') grid on; hold on; plot(x,exp(-x),'r:*');

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Other plotting commands: semilogx, semilogy

x = 0:0 1:10; 0:0.1:10; semilogx(10.^x,x)
10 9 8 7 6 5 4 3 2 1 0 0 10 10
0

x = 0:0 1:10; 0:0.1:10; semilogy(x, 10.^x)

10
10

10

10

10

10

10

10

10

10

10

10

10

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Other plotting commands: loglog

x = logspace(-1,2); loglog(x,exp(x),'-s') l l ( ( ) ) grid on
10
50

10

40

10

30

10

20

10

10

10

10

-1

10

10

10

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Matrices
A 2-D array, or matrix, of data is entered row by row, with spaces (or commas) separating entries within the row and semicolons separating the rows:
>> A = [1 2 3; 4 5 6; 7 8 9] A = 1 2 3 4 5 6 7 8 9

Extracting bits of matrices:

A(j,k) gives jth row, kth column A(2:3,1:2) A(2:3 1:2) rows 2 through 3 and columns 1 through 2 A([1,3], :) all of rows 1 and 3 A(:, 1) first column
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Numerical Array Concatenation

a=[1 2;3 4] a = 1 3 2 4 Use square brackets [ ] 4*a; 5*a, 6*a] 4 8 8 16 12 24

cat_a=[a, 2*a; 3*a, cat_a = 1 2 2 3 4 6 3 6 4 9 12 12 5 10 6 15 20 18

4*a

Use [ ] to combine existing arrays as matrix elements

CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

More on Matrices
zeros(n) zeros(m,n) a d( , ) rand(m,n) eye(m,n) ones(n) ones(m,n) size(A) length(A) Returns a n n matrix of zeros Returns a m n matrix of zeros Returns a m n matrix of random numbers Returns a m n Identity matrix Returns a n n matrix of ones Returns a m n matrix of ones For a m n matrix A, returns the row vector [m,n] containing the number of rows and columns in matrix d l i ti Returns the larger of the number of rows or columns in A
Dr. Tanvir Ahmed

CE 206: Engg. Computation Sessional

Diagonal Matrix
First define a vector containing the values of the diagonal entries (in order) then use the diag function to form the matrix
>> d = [ 3 4 2] D = di (d) [-3 2], diag(d) d = -3 4 2 D = -3 0 0 0 4 0 0 0 2

This command is useful when we need to construct very large matrices. If A is a matrix, diag(A) extracts the diagonal elements of the matrix.
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Sparse Matrix
These are generally large matrices that have a very small proportion of non-zero entries
>> i = [1, 3, 5]; j = [2,3,4]; >> v = [10 11 12]; >> S = sparse (i,j,v) S = ( , ) (1,2) 10 Creating a 5-by-4 sparse matrix S having (3,3) 11 only 3 non-zero values: (5,4) 12 S(1,2) = 10, >> T = full(S) S(3,3) = 11 and T = S(5,4) = 12 0 10 0 0 0 0 0 0 0 0 11 0 0 0 0 0 The full command displays the sparse 0 0 0 12 i matrix
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Matrix operations
Transpose Addition and Subtraction Scalar Multiplication Matrix Multiplication Matrix Inverse Matrix powers Determinant B = A C = A+B C = A-B

B = A where is a scalar A, C = A * B B = inv(A), A must be a square matrix B = A * A , A must be a square matrix det(A), A must be a square matrix
-standard standard -element-wise

Operators * ^ and / have two modes of operation

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Standard matrix product operation

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Elementwise matrix operations

To do element-wise operations, use the dot .(.*, ./, .^) BOTH dimensions must match (unless one is scalar) ).

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

for loop
Used when we want to repeat a segment of the code several times
Example: Test the assertion that the p ratio of the two successive values in the Fibonacci series approach the golden ratio of (5-1)/2 f . i.e. fn/fn-1 = (5-1)/2 where n = 3,4,5 ..

>> F(1) = 0; F(2) = 1; >> for i = 3:20 F(i) = F(i-1) + F(i-2); end >> plot(1:19, F(1:19)./F(2:20),'o' ) >> h ld on, xlabel('n') hold l b l(' ') >> plot(1:19, F(1:19)./F(2:20),'-' ) >> legend('Ratio of terms f_{n-1}/f_n') >> plot([0 20], (sqrt(5) 1)/2*[1 1] ' ') 20] (sqrt(5)-1)/2*[1,1],'--')
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

while loop
when we want to repeat a segment of the code several times until some logical condition is satisfied. The while loop is used when we do not know for certain how many iterations may be needed
Example: What is the greatest value of n that can be used in the sum and get a value less than 100? >> S = 1; n = 1; >> while S+ (n+1)^2 < 100 n = n+1; S = S + n^2; end d >> [n, S] ans = 6 91
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Logical and relational operators

== ~= > < >= <= & | && || Equal Not equal Greater than Less than Greater or equal G l Less or equal And Or x = -2.0000 3.1416 5.0000 -5.0000 -3.0000 -1.0000 >> x > 3 & x < 4 ans = 0 1 0 0 0 0 >> x > 3 | x == -3 ans = 0 1 1 0 1 0

Matlab represents true and false by means of the integers 0 and 1

CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Flow control using if/else/elseif

count=0; for n=1:length(x) if x(n)>0 count=count+1; end d end

CE 206: Engg. Computation Sessional

Example: Given x= sin(linspace(0 10*pi 100)) how many sin(linspace(0,10 pi,100)), of the entries are positive?

Dr. Tanvir Ahmed

The find command for vectors

returns a list of the positions (indices) of the elements of a vector satisfying a given condition.
Example: Produce a plot for y = e-x2 sin(3x) and mark all the points that have a value of y greater than 0.2 >> x = -1: 05:1; -1:.05:1; >> y = sin(3*pi*x).*exp(-x.^2); 1 >> k = find(y > 0.2) 0.8 k = 0.6 Columns 1 through 12 0.4 9 10 11 12 13 22 23 0.2 24 25 26 27 36 0 Columns 13 through 15 -0.2 37 38 39
-0.4 -0.6

>> plot(x,y,':'); hold on; -0.8 p ( ,y, ); ; 0.8 >> plot(x(k),y(k),'o) -1

CE 206: Engg. Computation Sessional

-1

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Dr. Tanvir Ahmed

Avoiding loops: efficient coding

Example 1: Given x= sin(linspace(0,10*pi,100)), how many of the entries are positive? count=0; for n=1:length(x) if x(n)>0 count=count+1; end end d

Using the find command find

Count=length(find(x>0));

Example 2 find the sum: 12 + 22 + 32 + . +1002 sum_sq=0; for n=1:100 sum_sq = sum_sq + n^2; end
CE 206: Engg. Computation Sessional

Sum_sq=sum((1:100).^2)

Vectorization
Dr. Tanvir Ahmed

Avoiding loops: efficient coding

Built-in functions (e.g. find, sum) will make it faster to write and execute
count=0; for n=1:length(x) g ( ) if x(n)>0 count=count+1; end end

Vectorized code is more efficient for MATLAB Use indexing and matrix operations to avoid loops

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Advanced graphics: plotting surfaces

- The command meshgrid is used to construct the (x, y) gridpoints h d d h ( ) d at certain intervals - Evaluate the function z = f(x y) at all the gridpoints f(x,y) -Use a surface plot feature (e.g. mesh,surf) to plot the surface
Example: Plot the surface defined by z = (x - 3)2 + (y 2)2 for 2 x 4 amd 1 y 3
Saddle

[X,Y] = meshgrid(2:.2:4, 1:.2:3); Z = (X-3).^2-(Y-2).^2; 1 mesh(X,Y,Z) 0.5 title('Saddle'), xlabel('x'),ylabel('y') 0

-0.5

-1 3 2.5 2 1.5 y 1 2 2.5 3.5 35 3 x 4

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Advanced graphics: surf and shading

[X,Y] [X Y] = meshgrid(2: 2:4 1: 2:3); meshgrid(2:.2:4, 1:.2:3); Z = (X-3).^2-(Y-2).^2; surf(X,Y,Z) title('Saddle) i l ( ddl ) shading faceted
1 0.5 0 -0.5 Saddle

shading flat colormap(gray) p(g y)

Saddle

-1 3 2.5 2 1.5 y 1 2 2.5 x 3.5 3 4

Saddle

0.5

0.5

-0.5

-0.5

-1 3 2.5 2 1.5 y 1 2 3 2.5 x 3.5 35 4

-1 3 2.5 2 1.5 y 1 2 2.5 x 3.5 35 3 4

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Advanced graphics: contour plot

Takes the same arguments as s rf or mesh surf
Example: Plot the contour of the surface defined by z = (sinx) (cosy) for x and - y x=-pi:0.1:pi; y=-pi:0.1:pi; [X,Y]=meshgrid(x,y); [X Y] h id( ) Z =sin(X).*cos(Y); contour(X,Y,Z,'LineWidth',2)
3 2 1 0 -1 -2

-3 -3 -2 -1 0 1 2 3

0.5

-0.5

-1 4 2 0 -2 -2 -4 -4 2 0 4

mesh(X,Y,Z); hold on; contour(X,Y,Z)

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Specialized plotting functions

polar-to make polar plots
polar(0:0.01:2 pi,cos((0:0.01:2 pi) 2)) polar(0:0.01:2*pi,cos((0:0.01:2*pi)*2))

bar-to make bar graphs

bar(1:10,rand(1,10)); bar(1:10 rand(1 10));

quiver-to add velocity vectors to a plot

[X,Y]=meshgrid(1:10,1:10); [X Y]=meshgrid(1:10 1:10); quiver(X,Y,rand(10),rand(10));

stairs plot stairs-plot piecewise constant functions

stairs(1:10,rand(1,10));

fill draws fill-draws and fills a polygon with specified vertices

fill([0 1 0.5],[0 0 1],'r');

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Scripting
Script files are ASCII (text) files containing MATLAB commands.

Commonly known as m-files because they have a .m extension

CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Function m-files
Example: Write a function m-file which calculates the area of a triangle having sides a, b and c using the formula: Area = s ( s a)( s b)( s c) Where s = (a+b+c)/2 The function name. Also the name of the m-file where the function definition will be stored. List of inputs Purpose of the function and how it can be used. Mainly y to aid the future users The code

List of output(s) p ( )

function [A] = area(a,b,c) % C Compute th area of a t i t the f triangle whose l h % sides have length a, b and c. % Inputs: % a,b,c: Lengths of sides % Output: % A: area of triangle g s = (a+b+c)/2; A = sqrt(s*(s-a)*(s-b)*(s-c)); %%%%%%%%% end of area %%%%%%%%%%%
CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Function overloading
MATLAB functions are generally overloaded Can take a variable number of inputs Can return a variable number of outputs You can overload your own functions by having variable input p g and output arguments
The following function plots a sine wave with frequency f1 on the range [0, 2], uses f1 as the input but displays a message when two inputs are given function plotSin(f1,f2) x linspace(0,2 pi,f1 16+1); x=linspace(0,2*pi,f1*16+1); figure if nargin == 1 plot(x,sin(f1*x)); plot(x sin(f1*x)); elseif nargin == 2 disp('Two inputs were given'); d end
Dr. Tanvir Ahmed

built in function nargin contains the number of i t inputs

CE 206: Engg. Computation Sessional