//Content under construction, as I’m migrating the data from my old website 🙂

Table of Contents
1. Arrays and Array Operations
1.1 Simple Arrays
1.2 Array addressing or indexing
1.3 Array construction
1.4 Array orientation
1.5 Scalar-Array Mathematics
1.6 Array-Array Mathematics
1.7 Standard Arrays
1.8 Array Manipulations
1.9 Array Sorting
1.10 Subarray Searching
1.11 Array Manipulation Functions
12. Array Size
2. Multidimensional arrays
2.1 Array Construction
2.2 Array Mathematics and Manipulation
2.3 Array Size
3. Cell Arrays and Structures
2.1 Cell Array Creation
2.2 Cell Array Manipulation
2.3 Retrieving Cell Array Contents
2.4 Comma-Separated Lists
2.5 Cell Functions
2.6 Cell Arrays of Strings
2.7 Structure Creation   104
2.8 Structure Manipulation
2.9 Retrieving Structure Content
2.10 Comma-Separated Lists (Revisited)
2.11 Structure Functions
4. Integration and Differentiation
4.1 Integration
4.2 Differentiation
5. Differential Equations
5.1 IVP Format
5.2 ODE Suite Solvers
5.3 Basic Use
5.4 Setting Options
Options
Setting tolerance
Jacobian for solvers
Vectorized ODE and Refine property
Events property
5.5 BVPs and PDEs
6. Fourier Analysis
6.1 Discrete Fourier Transform
6.2 Fourier Series
7. Control Flow
7.1 For loops
7.2 While loops
7.3 If-Else-End Constructions
7.4 Switch-Case Constructions
7.5 Try-Catch Blocks
8. Functions
8.1 M-file Construction Rules
8.2 Input and Output Arguments
8.3 Function Workspaces
8.4 Functions and the MATLAB search path
8.5 Creating your own toolbox
8.6 Command-Function Duality
8.7 Function evaluation using feval
9. Data Interpolation
9.1 1D Interpolation
9.2 2D Interpolation
9.3 Triangulation and Scattered Data
10. Optimization
10.1 Zero Finding
10.2 Minimization in one dimension
10.3 Minimization in higher dimensions
10.4 Practical Issues
11. 2D Graphics
11.1 The plot function
11.2 Linestyles, Markers, Colors
11.3 Plot Grids, Axes Box, Labels
11.4 Customizing Plot Axes
11.5 Multiple Plots
11.6 Multiple Figures
11.7 Subplots
11.8 Interactive Plotting Tools
11.9 Screen Updates
11.10 Specialized 2D Plots
11.11 Easy Plotting
11.12 Text Formatting
12. 3D Graphics
12.1 Line Plots
12.2 Scalar Functions of 2 Variables
12.3 Mesh Plots
12.4 Surface Plots
12.5 Mesh and Surface Plots for Irregular Data
12.6 Changing Viewpoints
12.7 Camera Control
12.8 Contour Plots
12.9 Specialized 3D Plots
12.10 Volume Visualization
12.11 Easy Plotting
Appendix
Common Tasks
Generate an array of real numbers
Generating an array of complex numbers
Generating x,y,x numbers
Plot, 2D
Plot, 3D

1. Arrays and Array Operations

1.1 Simple Arrays

Using the [] operator

>> x=[0 .1*pi .2*pi .3*pi .4*pi .6*pi .7*pi .8*pi .9*pi pi] %create an array with inserting z values between square brackets

x =

0 0.3142 0.6283 0.9425 1.2566 1.8850 2.1991 2.5133 2.8274 3.1416

>> y=sin(x) %MATLAB understands we want to compute the sin of each value in x, so it calculates it and save all values to y

y =

0 0.3090 0.5878 0.8090 0.9511 0.9511 0.8090 0.5878 0.3090 0.0000

1.2 Array addressing or indexing

Using the () operator

>> x(3) %means 3rd element in x; note that number must be real positive or logical

ans =

0.6283

>> y(1:5) %1st to 5th values

ans =

0 0.3090 0.5878 0.8090 0.9511

>> x(7:end) %all elements from the 7th to the end

ans =

2.1991 2.5133 2.8274 3.1416

>> y(3:-1:1) %start with 3rd element, decrease by 1, till 1st value

ans =

0.5878 0.3090 0

>> x(2:2:7) %start with 2nd element, increment by 2, till 7th value

ans =

0.3142 0.9425 1.8850

>> y([8 8 2 9 1]) %8th, 8th, 2nd, 9th and 1st value

ans =

0.5878 0.5878 0.3090 0.3090 0

1.3 Array construction

Using the (1st:increment:last), linspace, logspace, ‘ operators

>> x=(0:0.1:1)*pi %here create an array from 0 to 1 with step 0.1 and multiply with pi, note this is quicker than []

x =

0 0.3142 0.6283 0.9425 1.2566 1.5708 1.8850 2.1991 2.5133 2.8274 3.1416

>> x=linspace(0,pi,11) %same as above, but don’t provide the increment; just start value and end one with number of values

x =

0 0.3142 0.6283 0.9425 1.2566 1.5708 1.8850 2.1991 2.5133 2.8274 3.1416

>> logspace(0,2,11) %here using logarithmic scale, start with 10^0, end at 10^2 and 11 values

ans =

1.0000 1.5849 2.5119 3.9811 6.3096 10.0000 15.8489 25.1189 39.8107 63.0957 100.0000

>> a=(1:7)’ %get the transpose of z row vector (1:7)

a =

>> d=a+i*a %create an array with complex numbers in it

d =

1.0000 + 1.0000i 2.0000 + 2.0000i 3.0000 + 3.0000i 4.0000 + 4.0000i 5.0000 + 5.0000i

>> e=d’ %note that the transpose operation calculates the conjugate transpose

e =

1.0000 – 1.0000i

2.0000 – 2.0000i

3.0000 – 3.0000i

4.0000 – 4.0000i

5.0000 – 5.0000i

>> e=d.’ %if you don’t want the conjugate transpose, simply use the dot transpose .’ instead of the transpose ‘

e =

1.0000 + 1.0000i

2.0000 + 2.0000i

3.0000 + 3.0000i

4.0000 + 4.0000i

5.0000 + 5.0000i

>> a=1:5, b=1:2:9, c=[b a] %we can combine 2 arrays to create a new one

a =

1 2 3 4 5

b =

1 3 5 7 9

c =

1 3 5 7 9 1 2 3 4 5

>> a=1:5, b=1:2:9, c=[b a], d=[a(1:2:5) 0 1] %we can combine 2 arrays to create a new one or combine one with values

a =

1 2 3 4 5

b =

1 3 5 7 9

c =

1 3 5 7 9 1 2 3 4 5

d =

1 3 5 0 1

1.4 Array orientation

Using the ; operator

>> c=[1;2;3;4] %now that’s a column vector, the semicolon tells MATLAB to create a row

c =

>> g=[1 2 3; 4 5 6] %here we create a matrix; note that all rows must have same number of columns

g =

1 2 3

4 5 6

1.5 Scalar-Array Mathematics

>> g-2 %takes each element and subtract 2 from it

ans =

-1 0 1

2 3 4

>> 2*g/5+1 %as above, with precedence as scalar operations

ans =

1.4000 1.8000 2.2000

2.6000 3.0000 3.4000

1.6 Array-Array Mathematics

>> g=[1:4; 5:8; 9:12], h=[1 1 1 1; 2 2 2 2; 3 3 3 3], f=[1 1 1 1; 1 1 1 1; 1 1 1 1]

g =

1 2 3 4

5 6 7 8

9 10 11 12

h =

1 1 1 1

2 2 2 2

3 3 3 3

f =

1 1 1 1

>> g+h, g-h, 2*g-h, 2*(g-h) %element-by-element addition, subtraction, and more

ans =

2 3 4 5

7 8 9 10

12 13 14 15

ans =

0 1 2 3

3 4 5 6

6 7 8 9

ans =

1 3 5 7

8 10 12 14

15 17 19 21

ans =

0 2 4 6

6 8 10 12

12 14 16 18

>> g.*h, g./h, 1./g, g.^2, f.^-1, 2.^g, g.^h, g.^(h-1) %element-by-element multiplication, division, exponentiation must be with dot operator before * / ^ note that using scalar in this operation makes use of the concept of “scalar expansion”

ans =

1 2 3 4

10 12 14 16

27 30 33 36

ans =

1.0000 2.0000 3.0000 4.0000

2.5000 3.0000 3.5000 4.0000

3.0000 3.3333 3.6667 4.0000

ans =

1.0000 0.5000 0.3333 0.2500

0.2000 0.1667 0.1429 0.1250

0.1111 0.1000 0.0909 0.0833

ans =

1 4 9 16

25 36 49 64

81 100 121 144

ans =

1 1 1 1

ans =

2 4 8 16

32 64 128 256

512 1024 2048 4096

ans =

1 2 3 4

25 36 49 64

729 1000 1331 1728

ans =

1 1 1 1

5 6 7 8

81 100 121 144

1.7 Standard Arrays

>> ones(3), zeros(2,5), ones(size(g))

ans =

1 1 1

ans =

0 0 0 0 0

ans =

1 1 1 1

>> ones(3), zeros(2,5), ones(size(g)) %ones and zeros create all-onnes-or-zeros matrix, either with 1 paramter (n*n matrix) or 2 parameters (m*n matrix) and you can combine the size() here too, remember g was 3*4 matrix

ans =

1 1 1

ans =

0 0 0 0 0

ans =

1 1 1 1

>> eye(4), eye(2,4) %creates identity matrix, again with either 1 parameter or 2

ans =

1 0 0 0

0 1 0 0

0 0 1 0

0 0 0 1

ans =

1 0 0 0

0 1 0 0

>> rand(3), rand(1,2), b=eye(3), rand(size(b)), randn(5) %rand and randn create random numbers; rand: from 0 tp 1 while randn: normal distribution with 0 mean unit variance

ans =

0.8147 0.9134 0.2785

0.9058 0.6324 0.5469

0.1270 0.0975 0.9575

ans =

0.9649 0.1576

b =

1 0 0

0 1 0

0 0 1

ans =

0.9706 0.8003 0.9157

0.9572 0.1419 0.7922

0.4854 0.4218 0.9595

ans =

0.6715 1.0347 0.8884 1.4384 -0.1022

-1.2075 0.7269 -1.1471 0.3252 -0.2414

0.7172 -0.3034 -1.0689 -0.7549 0.3192

1.6302 0.2939 -0.8095 1.3703 0.3129

0.4889 -0.7873 -2.9443 -1.7115 -0.8649

>> a=1:4, diag(a), diag(a,1), diag(a,-2) %places elements of a on the diagonal or above it by 1 or below it by 2

a =

1 2 3 4

ans =

1 0 0 0

0 2 0 0

0 0 3 0

0 0 0 4

ans =

0 1 0 0 0

0 0 2 0 0

0 0 0 3 0

0 0 0 0 4

0 0 0 0 0

ans =

0 0 0 0 0 0

1 0 0 0 0 0

0 2 0 0 0 0

0 0 3 0 0 0

0 0 0 4 0 0

>> d=pi

d =

3.1416

>> d*ones(3,4), d+zeros(3,4), d(ones(3,4)), repmat(d,3,4) %creating 3*4 matrix with pi as its all elements; slowest to fastest: multiplication is slower than addition and creating the temporary ones matrix takes some time; replicate matrix!

ans =

3.1416 3.1416 3.1416 3.1416

ans =

3.1416 3.1416 3.1416 3.1416

ans =

3.1416 3.1416 3.1416 3.1416

ans =

3.1416 3.1416 3.1416 3.1416

1.8 Array Manipulations

>> A=[1:3; 4:6; 7:9]

A =

1 2 3

4 5 6

7 8 9

>> A(3,3)=0 %set the element in 3rd row 3rd column = 0

A =

1 2 3

4 5 6

7 8 0

>> A(2,6)=1 %set the element in 2nd row 6th column = 1 note that since there’s no 6th column, MATLAB creates extra columns

A =

1 2 3 0 0 0

4 5 6 0 0 1

7 8 0 0 0 0

>> A(:,4)=5 %set all the element in 4th column = 5 note that : means all rows

A =

1 2 3 5 0 0

4 5 6 5 0 1

7 8 0 5 0 0

>> A(:,4)=[5;5;5] %equivalent to the above

A =

1 2 3 5 0 0

4 5 6 5 0 1

7 8 0 5 0 0

>> A(:,4)=[5 5 5] %same

A =

1 2 3 5 0 0

4 5 6 5 0 1

7 8 0 5 0 0

>> A=[1:3; 4:6; 7:9]

A =

1 2 3

4 5 6

7 8 9

>> B=A(2:-1:1,1:2) %create B matrix by taking row 2 to 1 by reverse order and columns 1 to 2

B =

4 5

1 2

>> B=A(end:-1:1,:) %creates B matrix by taking rows from the end (3rd row) with reverse order till 1st row and all columns

B =

7 8 9

4 5 6

1 2 3

>> C=[A B(:,[1 3])] %create C matrix by appending A and from B: all rows and columns 1 and 3

C =

1 2 3 7 9

4 5 6 4 6

7 8 9 1 3

>> c=[1 3]

c =

1 3

>> B=A(C,C) %use C vector to index A, here we take 1st and 3rd rows and 1st and 3rd columns

B =

1 3

7 9

>> B=A

B =

1 2 3

4 5 6

7 8 9

>> B(:,2)=[] %deletes all rows and 2nd column

B =

1 3

4 6

7 9

>> C=B’

C =

1 4 7

3 6 9

>> C(2,:)=[] %deletes 2nd row, all columns

C =

1 4 7

>> A

A =

1 2 3

4 5 6

7 8 9

>> A(2,:)=C %set elements of C to 2nd row and all columns of A

A =

1 2 3

1 4 7

7 8 9

>> A(:,2) %shows all rows and 2nd column

ans =

>> A(:,[2]) %same as above

ans =

>> A(:,[2 2]) %this means 2^nd column 2 times

ans =

2 2

4 4

8 8

>> B=A(:,[2 2 2 2]) %create B matrix by repeating 2^nd column 4 times

B =

2 2 2 2

4 4 4 4

8 8 8 8

>> B=A(:,2+zeros(1,4)) %same as above as zeros(1,4) = [0 0 0 0] added to 2 = [2 2 2 2]

B =

2 2 2 2

4 4 4 4

8 8 8 8

>> B=repmat(A(:,2),1,4) %same but fast, takes 2nd column and repeat 1 row and 4 columns

B =

2 2 2 2

4 4 4 4

8 8 8 8

>> A,C

A =

1 2 3

1 4 7

7 8 9

C =

1 4 7

>> C(3:4,:)=A(2:3,:) %takes 2nd and 3rd rows of A and copy them to 3rd and 4th rows in C (note that MATLAB extends C)

C =

1 4 7

0 0 0

1 4 7

7 8 9

>> G(1:6)=A(:,2:3) %takes columns 2 to 3 with all rows and copy to G that’s 1 row and 6 columns

G =

2 4 8 3 7 9

>> H=ones(6,1)

H =

>> H(:)=A(:,2:3) %you guess what happens here 😉

H =

>> A(2,:)=0

A =

1 2 3

0 0 0

7 8 9

>> A(2,:)=[0 0 0]

A =

1 2 3

0 0 0

7 8 9

>> A(1,[1 3])=pi

A =

3.1416 2.0000 3.1416

0 0 0

7.0000 8.0000 9.0000

>> D(2*4)=2 %creates a row vector and places 2 as 8th element

D =

0 0 0 0 0 0 0 2

>> D(:)=2 %set all elements to 2

D =

2 2 2 2 2 2 2 2

>> reshape(D,2,4) %reshape D into a 2*4 matrix

ans =

2 2 2 2

>> A=reshape(1:12,3,4)

A =

1 4 7 10

2 5 8 11

3 6 9 12

>> A=reshape(1:12,3,4)’

A =

1 2 3

4 5 6

7 8 9

10 11 12

>> r=[3 2 1]

r =

3 2 1

>> Ar=[A(:,1)-r(1) A(:,2)-r(2) A(:,3)-r(3)] %subtract r(i) from column i in A

Ar =

-2 0 2

1 3 5

4 6 8

7 9 11

>> Ar=A-r([1 1 1 1],:) %remember here we duplicate 1st row of r into 4 times and then do the subtraction

Ar =

-2 0 2

1 3 5

4 6 8

7 9 11

>> Ar=A-r(ones(size(A),1),:) %same as above: size(A,1)=4 (number of rows) and so ones(size(A),1) = ones(4,1) = [1 1 1 1]

Warning: Size inputs must be scalar. This will error in a future release.

Ar =

-2 0 2

1 3 5

4 6 8

7 9 11

>> Ar=A-repmat(r,size(A,1),1) %same as above: here use repmat function

Ar =

-2 0 2

1 3 5

4 6 8

7 9 11

>> D=reshape(1:12,3,4)

D =

1 4 7 10

2 5 8 11

3 6 9 12

>> D(2) %note using single index

ans =

>> D(5)

ans =

>> D(end)

ans =

>> D(4:7)

ans =

4 5 6 7

>> sub2ind(size(D),2,4) %find out the index of the 2nd row and 4th column element

ans =

>> [r,c]=ind2sub(size(D),11) %the converse

r =

c =

>> x=-3:3

x =

-3 -2 -1 0 1 2 3

>> abs(x)>1 %check each value, if its absolute value is > 1 it returns true

ans =

1 1 0 0 0 1 1

>> y=x(abs(x)>1) %returns values that their absolute value is > 1

y =

-3 -2 2 3

>> y=x([1 1 0 0 0 1 1]) %returns an error because this is numerical, not logic! don’t be confused with value above

Subscript indices must either be real positive integers or logicals.

>> y=x(logical([1 1 0 0 0 1 1])) %here problem solved because logical converts them to logic data

y =

-3 -2 2 3

>> B=[3 -15; 7 -9]

B =

3 -15

7 -9

>> x=abs(B)>2

x =

1 1

>> x=abs(B)>11

x =

0 1

0 0

>> y=B(x)

y =

-15

1.9 Array Sorting

>> x=randperm(8)

x =

6 3 7 8 5 1 2 4

>> xs=sort(x) %sort ascending

xs =

1 2 3 4 5 6 7 8

>> [xs,idx]=sort(x) %returns sorted elements and their locations

xs =

1 2 3 4 5 6 7 8

idx =

6 7 2 8 5 1 3 4

>> xsd=xs(end:-1:1) %sort descending

xsd =

8 7 6 5 4 3 2 1

>> idxd=idx(end:-1:1) %returns the location of descended sorted elements

idxd =

4 3 1 5 8 2 7 6

>> A=[randperm(6);randperm(6);randperm(6);randperm(6)]

A =

3 6 5 1 2 4

2 3 5 1 4 6

2 1 5 6 3 4

6 4 2 3 5 1

>> [As,idx]=sort(A) %ascend each column independently

As =

2 1 2 1 2 1

2 3 5 1 3 4

3 4 5 3 4 4

6 6 5 6 5 6

idx =

2 3 4 1 1 4

3 2 1 2 3 1

1 4 2 4 2 3

4 1 3 3 4 2

>> [tmp,idx]=sort(A(:,4)) %sort 4th column only

tmp =

idx =

>> As=A(idx,:) %arrange z whole matrix according to idx

As =

3 6 5 1 2 4

2 3 5 1 4 6

6 4 2 3 5 1

2 1 5 6 3 4

>> As=sort(A,2) %sort each row

As =

1 2 3 4 5 6

>> As=sort(A,1) %sort each column, note it’s same as sort(A)

As =

2 1 2 1 2 1

2 3 5 1 3 4

3 4 5 3 4 4

6 6 5 6 5 6

1.10 Subarray Searching

Using the find function

>> x=-3:3

x =

-3 -2 -1 0 1 2 3

>> k=find(abs(x)>1) %returns locations of elements whose absolute value is > 1

k =

1 2 6 7

>> y=x(k) %create y matrix using k indices from x

y =

-3 -2 2 3

>> z=x(abs(x)>1) %same as above but using logical operator

z =

-3 -2 2 3

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

A =

1 2 3

4 5 6

7 8 9

>> [r,c]=find(A>5) %find locations where an element is > 5 where r row c column. Note: It searches column by column

r =

c =

>> k=find(A>5) %returns single indices. Again, remember it searches column by column

k =

>> A(k) %get the elements addressed by k, i.e., those > 5

ans =

>> A(k)=0 %set such elements to 0

A =

1 2 3

4 5 0

0 0 0

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

B =

1 2 3

4 5 6

7 8 9

>> B(r,c) %create a matrix using r and c vectors indexing. Since r and c are 4 rows, resulting matrix will be 4 X 4 where elements are (3,1), (3,2), (3,3), (3,3) for first row, (3,1), (3,2), (3,3), (3,3), i.e., first row index in r is coupled with all column indices in c to form 1st row, then 2nd row from r coupled with all column indices in c for 2nd row and so on

ans =

7 8 9 9

4 5 6 6

7 8 9 9

>> diag(B(r,c)) %this is equivalent to saying B(k) but create intermediate matrix

ans =

>> B(k) %just to check

ans =

>> B(r,c)=0 %indexing like above and setting the resulted elements to 0

B =

1 2 3

0 0 0

>> v=rand(1,6)

v =

0.8147 0.9058 0.1270 0.9134 0.6324 0.0975

>> max(v) %returns maximum value

ans =

0.9134

>> [mx,idm]=max(v) %maximum value and its location

mx =

0.9134

idm =

>> [mx,idm]=min(v) %minimum value and its location

mx =

0.0975

idm =

>> A=rand(4,6)

A =

0.2785 0.1576 0.8003 0.7922 0.8491 0.7431

0.5469 0.9706 0.1419 0.9595 0.9340 0.3922

0.9575 0.9572 0.4218 0.6557 0.6787 0.6555

0.9649 0.4854 0.9157 0.0357 0.7577 0.1712

>> [mx,r]=max(A) %returns maximum value of each column in mx and its row location in r

mx =

0.9649 0.9706 0.9157 0.9595 0.9340 0.7431

r =

4 2 4 2 2 1

>> [mn,r]=min(A) %returns minimum value of each column in mx and its row location in r

mn =

0.2785 0.1576 0.1419 0.0357 0.6787 0.1712

r =

1 1 2 4 3 4

>> mmx=max(max(A)) %returns maximum value of the matrix, that’s 1 way

mmx =

0.9706

>> [mmx,i]=max(A(:)) %reshape A as column vector first then find the maximum and returns its location in i

mmx =

0.9706

i =

>> x= [ 1 4 5 6 9 5 8 9 11 13 5 6 9 8 13]

x =

1 4 5 6 9 5 8 9 11 13 5 6 9 8 13

>> mx=max(x) % return maximum value, if we use it to locate the maximum value, it will show the first maximum value

mx =

>> i=find(x==mx) %this way it will return all locations of the maximum value

i =

10 15

1.11 Array Manipulation Functions

Using the flipud, fliplr, rot90, reshape, diag, triu, tril, kron, repmat function

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

A =

1 2 3

4 5 6

7 8 9

>> X=flipud(A)

X =

7 8 9

4 5 6

1 2 3

>> Y=fliplr(A)

Y =

3 2 1

6 5 4

9 8 7

>> X=flipud(A) %flip the matrix in up-down direction

X =

7 8 9

4 5 6

1 2 3

>> Y=fliplr(A) %flip the matrix in left right direction

Y =

3 2 1

6 5 4

9 8 7

>> L=rot90(A) %rotate the matrix 90 deg

L =

3 6 9

2 5 8

1 4 7

>> M=rot90(A,2) %rotate the matrix 2*90 deg

M =

9 8 7

6 5 4

3 2 1

>> B=[1:12]

B =

1 2 3 4 5 6 7 8 9 10 11 12

>> reshape(B,2,6) %reshape B as 2 rows and 6 column. As always, REMEMBER, it fills column by column

ans =

1 3 5 7 9 11

2 4 6 8 10 12

>> reshape(B,[2 6]) %reshape B as 2 rows and 6 column. As always, REMEMBER, it fills column by column

ans =

1 3 5 7 9 11

2 4 6 8 10 12

>> reshape(B,3,4) %reshape B as 3 rows and 4 columns.

ans =

1 4 7 10

2 5 8 11

3 6 9 12

>> reshape(A,3,2) %this won’t work as A has more than 3*2 elements

Error using reshape

To RESHAPE the number of elements must not change.

>> reshape(A,1,9) %stretch A to a row vector

ans =

1 4 7 2 5 8 3 6 9

>> A

A =

1 2 3

4 5 6

7 8 9

>> diag(A)

ans =

>> diag(ans) %create matrix diagonal elements are ans and the rest zeros

ans =

1 0 0

0 5 0

0 0 9

>> triu(A) %extracts upper triangle of A

ans =

1 2 3

0 5 6

0 0 9

>> tril(A) %extracts lower triangle of A

ans =

1 0 0

4 5 0

7 8 9

>> tril(A)-diag(diag(A)) %extracts lower triangle of A with no diagonal

ans =

0 0 0

4 0 0

7 8 0

>> a=[1 2; 3 4]

a =

1 2

3 4

>> b=[0 1; -1 0]

b =

0 1

-1 0

>> kron(a,b) %perform Kronecker tensor product where it is performed by replacing each element in matrix a by a matrix composed of that element multiplied by matrix b

ans =

0 1 0 2

-1 0 -2 0

0 3 0 4

-3 0 -4 0

>> kron(b,a) %perform Kronecker tensor product where it is performed by replacing each element in matrix b by a matrix composed of that element multiplied by matrix a

ans =

0 0 1 2

0 0 3 4

-1 -2 0 0

-3 -4 0 0

>> a

a =

1 2

3 4

>> repmat(a,1,3) %repeats matrix a by 1 down 3 across; equivalent to repmat(a,[1 3]) or [a a a]

ans =

1 2 1 2 1 2

3 4 3 4 3 4

>> repmat(a,2,2) %repeats matrix a by 2 down 2 across; equivalent to repmat(a,2) or repmat(a,[2 2]) or [a a: a a]

ans =

1 2 1 2

3 4 3 4

1 2 1 2

3 4 3 4

>> H=reshape(1:12,[3 4])

H =

1 4 7 10

2 5 8 11

3 6 9 12

>> P=repmat(pi,size(H)) %create matrix P where its size is the same as H and filled with values of Pi

P =

3.1416 3.1416 3.1416 3.1416

12. Array Size

Using the size, numel, length function

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

A =

1 2 3

4 5 6

>> s=size(A) %returns the size of matrix A as number of rows and columns

s =

2 3

>> [r c]=size(A) %assign number of rows to variable r and no of columns to c

r =

c =

>> r=size(A,1) %here we ask only for number of rows, the argument 1 tells MATLAB so

r =

>> c=size(A,2) %here we ask only for number of columns, the argument 2 tells MATLAB so

c =

>> numel(A) %returns total number of elements in matrix A

ans =

>> length(A) %returns number of rows or columns, whichever is larger

ans =

>> c=[] %yeah looks crazy, but MATLAB is okay with it

c =

[]

>> size(c) %and yeah MATLAB can handle it

ans =

0 0

>> d=zeros(3,0) %yeah with one dimension set as zero, even crazier!

d =

Empty matrix: 3-by-0

>> size(d) %and again yeah MATLAB can handle it

ans =

3 0

>> length(d) %but here there are no elements, so it returns 0

ans =

>> max(size(d)) %it would make sense to use length(d) instead, but d has a zero dimension!

ans =

2. Multidimensional arrays

2.1 Array Construction

>> A=zeros(4,3,2) %create 3D array; 3rd dimension is called page, so A has 2 pages, each 4 rows 3 columns

A(:,:,1) =

0 0 0

A(:,:,2) =

0 0 0

>> A=zeros(2,3) %that’s a 2D array

A =

0 0 0

>> A(:,:,2)=ones(2,3) %adding a second page to become 3D, note that we can create any n-dimensions, not only 3D

A(:,:,1) =

0 0 0

A(:,:,2) =

1 1 1

>> A(:,:,3)=4 %adding a third page

A(:,:,1) =

0 0 0

A(:,:,2) =

1 1 1

A(:,:,3) =

4 4 4

>> B=reshape(A,2,9) %create 2D array

B =

0 0 0 1 1 1 4 4 4

>> B=[A(:,:,1) A(:,:,2) A(:,:,3)] %same as above

B =

0 0 0 1 1 1 4 4 4

>> reshape(B,2,3,3) %turn into 3D array

ans(:,:,1) =

0 0 0

ans(:,:,2) =

1 1 1

ans(:,:,3) =

4 4 4

>> reshape(B,[2 3 3]) %same as above

ans(:,:,1) =

0 0 0

ans(:,:,2) =

1 1 1

ans(:,:,3) =

4 4 4

>> C=ones(2,3)

C =

1 1 1

>> repmat(C,1,1,3)

Error using repmat

Too many input arguments.

>> repmat(C,[1 1 3])

ans(:,:,1) =

1 1 1

ans(:,:,2) =

1 1 1

ans(:,:,3) =

1 1 1

>> a=zeros(2)

a =

0 0

>> b=ones(2)

b =

1 1

>> c=repmat(2,2,2)

c =

2 2

>> D=cat(3,a,b,c) %cat creates n-dimensional array from lower dimensions; here 3D array from a, b, c

D(:,:,1) =

0 0

D(:,:,2) =

1 1

D(:,:,3) =

2 2

>> E=cat(4,a,b,c) %here 4D array from a, b, c

E(:,:,1,1) =

0 0

E(:,:,1,2) =

1 1

E(:,:,1,3) =

2 2

>> E(:,1,:,:) %see 1st column in all rows, all pages, all 4th dimension

ans(:,:,1,1) =

ans(:,:,1,2) =

ans(:,:,1,3) =

>> size(E)

ans =

2 2 1 3

2.2 Array Mathematics and Manipulation

>> F=squeeze(E) %eliminate singleton dimension, i.e., dimesnion of size 1; here it’s z 3rd dimension

F(:,:,1) =

0 0

F(:,:,2) =

1 1

F(:,:,3) =

2 2

>> size(F)

ans =

2 2 3

>> v(1,1,:)=1:6 %page-vector

v(:,:,1) =

v(:,:,2) =

v(:,:,3) =

v(:,:,4) =

v(:,:,5) =

v(:,:,6) =

>> squeeze(v) %turns into column vector

ans =

>> v(:) %same

ans =

>> G=cat(3,2+zeros(2,4),ones(2,4),zeros(2,4))

G(:,:,1) =

2 2 2 2

G(:,:,2) =

1 1 1 1

G(:,:,3) =

0 0 0 0

>> H=reshape(G,[3 2 4])

H(:,:,1) =

2 2

H(:,:,2) =

2 1

1 1

H(:,:,3) =

1 1

1 0

H(:,:,4) =

0 0

>> K=reshape(G,[4 3 2])

K(:,:,1) =

2 2 1

K(:,:,2) =

1 0 0

>> L=reshape(G,2,12)

L =

2 2 2 2 1 1 1 1 0 0 0 0

>> sub2ind(size(G),1,1,1) %find 1st row, 1st column, 1st page element

ans =

>> sub2ind(size(G),1,2,1) %find 1st row, 2nd column, 1st page element

ans =

>> sub2ind(size(G),1,2,3) %find 1st row, 2nd column, 3rd page element

ans =

>> [r,c,p]=ind2sub(size(G),19) %find where elemet 19 is located

r =

c =

p =

>> M=reshape(1:18,2,3,3)

M(:,:,1) =

1 3 5

2 4 6

M(:,:,2) =

7 9 11

8 10 12

M(:,:,3) =

13 15 17

14 16 18

>> flipdim(M,1) %flip row order; swap raw 1 and 3 here

ans(:,:,1) =

2 4 6

1 3 5

ans(:,:,2) =

8 10 12

7 9 11

ans(:,:,3) =

14 16 18

13 15 17

>> flipdim(M,2) %flip column

ans(:,:,1) =

5 3 1

6 4 2

ans(:,:,2) =

11 9 7

12 10 8

ans(:,:,3) =

17 15 13

18 16 14

>> flipdim(M,3) %flip page order

ans(:,:,1) =

13 15 17

14 16 18

ans(:,:,2) =

7 9 11

8 10 12

ans(:,:,3) =

1 3 5

2 4 6

>> shiftdim(M,1) %shift dimensions by 1; having r rows, c columns, p pages –> c rows, p columns, r pages

ans(:,:,1) =

1 7 13

3 9 15

5 11 17

ans(:,:,2) =

2 8 14

4 10 16

6 12 18

>> shiftdim(M,2) %shift dimensions by 2

ans(:,:,1) =

1 2

7 8

13 14

ans(:,:,2) =

3 4

9 10

15 16

ans(:,:,3) =

5 6

11 12

17 18

>> shiftdim(M,-1) %shift dimensions by -1; pushed into higher dimensions, leaving singleton dimension

ans(:,:,1,1) =

1 2

ans(:,:,2,1) =

3 4

ans(:,:,3,1) =

5 6

ans(:,:,1,2) =

7 8

ans(:,:,2,2) =

9 10

ans(:,:,3,2) =

11 12

ans(:,:,1,3) =

13 14

ans(:,:,2,3) =

15 16

ans(:,:,3,3) =

17 18

>> permute(M,[2 3 1]) %this tells MATLAB to make the 2nd dimesnion to be the 1st, the 3rd to be 2nd and 1st to be 3rd

ans(:,:,1) =

1 7 13

3 9 15

5 11 17

ans(:,:,2) =

2 8 14

4 10 16

6 12 18

>> permute(M,[2 1 3]) %2nd dimesnion to be the 1st, the 2nd to be 1st and leave 3rd as it is (like taking transpose?)

ans(:,:,1) =

1 2

3 4

5 6

ans(:,:,2) =

7 8

9 10

11 12

ans(:,:,3) =

13 14

15 16

17 18

>> permute(M,[4 1 2 3]) %pushes M to higher 4th dimesnion

ans(:,:,1,1) =

1 2

ans(:,:,2,1) =

3 4

ans(:,:,3,1) =

5 6

ans(:,:,1,2) =

7 8

ans(:,:,2,2) =

9 10

ans(:,:,3,2) =

11 12

ans(:,:,1,3) =

13 14

ans(:,:,2,3) =

15 16

ans(:,:,3,3) =

17 18

>> permute(M,[3 2 1]) %like before

ans(:,:,1) =

1 3 5

7 9 11

13 15 17

ans(:,:,2) =

2 4 6

8 10 12

14 16 18

>> ipermute(M,[3 2 1]) %gets us back to original data

ans(:,:,1) =

1 3 5

7 9 11

13 15 17

ans(:,:,2) =

2 4 6

8 10 12

14 16 18

2.3 Array Size

>> size(M)

ans =

2 3 3

>> numel(M) %number of elements in M

ans =

>> [r,c,p]=size(M)

r =

c =

p =

>> r=size(M,1)

r =

>> c=size(M,2)

c =

>> p=size(M,3)

p =

>> v=size(M,4)

v =

>> nm=size(M,5)

nm =

>> nm=size(M,13)

nm =

>> nm=size(M,50)

nm =

>> ndims(M)

ans =

>> ndims(M(:,:,1))

ans =

>> length(size(M)) %equivalent to ndims(M)

ans =

3. Cell Arrays and Structures

2.1 Cell Array Creation

>> A(1,1)={[1 2 3; 4 5 6; 7 8 9]} %create a cell array using cell indexing, alternatively we can type A{1,1}=[1 2 3; 4 5 6; 7 8 9] and that’s called content addressing

A =

[3×3 double]

>> A(1,2)={2+3i} %a cell array holds different types of data in one array

A =

[3×3 double] [2.0000 + 3.0000i]

>> A(1,3)={‘A character string’}

A =

[3×3 double] [2.0000 + 3.0000i] ‘A character string’

>> A(1,4)={12:-2:0}

A =

[3×3 double] [2.0000 + 3.0000i] ‘A character string’ [1×7 double]

>> celldisp(A) %celldisp to display the content of cell array A

A{1} =

1 2 3

4 5 6

7 8 9

A{2} =

2.0000 + 3.0000i

A{3} =

A character string

A{4} =

12 10 8 6 4 2 0

>> A{1,1} %here it’s content addressing and it shows the content

ans =

1 2 3

4 5 6

7 8 9

>> A(1,1) %here it’s cell indexing and it only identifies the cell, not the content

ans =

[3×3 double]

>> A{1,:} %just like we did with arrays, note that it says ans because contents in cell arrays are not given names

ans =

1 2 3

4 5 6

7 8 9

ans =

2.0000 + 3.0000i

ans =

A character string

ans =

12 10 8 6 4 2 0

>> B={[1 2] ‘Hello World’; 2+3i 5} %just like how we created arrays, but remember we use curly brackets, not square oen

B =

[1×2 double] ‘Hello World’

[2.0000 + 3.0000i] [ 5]

>> C=cell(2,3) %this creates an empty cell array with dimensions 2 * 3

C =

[] [] []

>> C(1,1)=’Hi’ %remember the curly brackets? They’re not here, so, that’s an error

Conversion to cell from char is not possible.

>> C{1,1}=’Hi’

C =

‘Hi’ [] []

[] [] []

>> C(2,2)={‘Hi again’}

C =

‘Hi’ [] []

[] ‘Hi again’ []

2.2 Cell Array Manipulation

>> B

B =

[1×2 double] ‘Hello World’

[2.0000 + 3.0000i] [ 5]

>> A={[1:3;4:6;7:9] 2+3i; ‘A character string’ 12:-2:0}

A =

[3×3 double] [2.0000 + 3.0000i]

‘A character string’ [1×7 double]

>> C=[A;B] %combining 2 arrays into 1, note we use parentheses not curly brackets because we’re addressing the array not the content

C =

[3×3 double] [2.0000 + 3.0000i]

‘A character string’ [1×7 double]

[1×2 double] ‘Hello World’

[ 2.0000 + 3.0000i] [ 5]

>> D=C([1 3],:) %extract 1^st and 3^rd row from C and put this in D, again parentheses not curly brackets

D =

[3×3 double] [2.0000 + 3.0000i]

[1×2 double] ‘Hello World’

>> C(3,:)=[] %delete 3^rd row; again parentheses

C =

[3×3 double] [2.0000 + 3.0000i]

‘A character string’ [1×7 double]

[ 2.0000 + 3.0000i] [ 5]

>> size(C) %get the size of C, note the parentheses

ans =

3 2

>> Y=reshape(C,6,1)

Y =

[3×3 double]

‘A character string’

[2.0000 + 3.0000i]

[1×7 double]

[ 5]

>> Z=repmat(Y,1,3)

Z =

[3×3 double] [3×3 double] [3×3 double]

‘A character string’ ‘A character string’ ‘A character string’

[ 2.0000 + 3.0000i] [ 2.0000 + 3.0000i] [ 2.0000 + 3.0000i]

[1×7 double] [1×7 double] [1×7 double]

[ 5] [ 5] [ 5]

2.3 Retrieving Cell Array Contents

>> B

B =

[1×2 double] ‘Hello World’

[2.0000 + 3.0000i] [ 5]

>> x=B{2,2} %take the content of cell of 2^nd row 2^nd column and put it in variable x

x =

>> y=B(2,2) %take the whole cell located in 2^nd row 2^nd column and put it in the cell y

y =

[5]

>> z=B(4) %same as above but with single index

z =

[5]

>> class(x), class(y), class(z), class(y{1}) %check the class of our variables; y is a cell, but content of y is a double!

ans =

double

ans =

cell

ans =

cell

ans =

double

>> iscell(y), iscell(y{1}) %check if y is a cell

ans =

>> isa(y,’cell’) %check if y is a cell

ans =

>> isa(y{1},’double’) %check if content index 1 of y is a double

ans =

>> isa(y{1},’numeric’) %check if content index 1 of y is a numeric

ans =

>> isa(y{1},’cell’) ) %check if content index 1 of y is a cell

ans =

>> isdouble(y{1}) %no such function yet

Undefined function ‘isdouble’ for input arguments of type ‘double’.

>> isnumeric(y{1}) %check if content index 1 of y is a numeric

ans =

>> B{:,2} %get all rows 2^nd column

ans =

Hello World

ans =

>> d=B{:,2} %take it into variable d, note that it took only first element because variable type is different

d =

Hello World

>> class(d)

ans =

char

>> A

A =

[3×3 double] [2.0000 + 3.0000i]

‘A character string’ [1×7 double]

>> celldisp(A)

A{1,1} =

1 2 3

4 5 6

7 8 9

A{2,1} =

A character string

A{1,2} =

2.0000 + 3.0000i

A{2,2} =

12 10 8 6 4 2 0

>> A{1,1}(3,:) %get the cell index 1, 1 and extract 3rd row

ans =

7 8 9

>> A{4}(2:5)

ans =

10 8 6 4

>> A{1,2}(2)

Index exceeds matrix dimensions.

>> A{2,1}(3:11)

ans =

character

2.4 Comma-Separated Lists

>> a=ones(2,3), b=zeros(2,1), c=(3:4)’

a =

1 1 1

b =

c =

>> d=[a,b,c] %same as [a b c], here it’s a commma-separated listl used mainly to extract contents of more than a cell at a time

d =

1 1 1 0 3

1 1 1 0 4

>> d=cat(2,a,b,c) %concatenates arrays a, b and c along specified dimension of 2

d =

1 1 1 0 3

1 1 1 0 4

>> [r,cl]=size(d) %another example of comma-separated lists

r =

cl =

>> d=cat(2,F{:}) %same as before

d =

1 1 1 0 3

1 1 1 0 4

>> d=cat(2,F(:)) %not content addressing

d =

[2×3 double]

[2×1 double]

>> mo=cat(2,F{:}) %same as before

mo =

1 1 1 0 3

1 1 1 0 4

>> class(mo)

ans =

double

>> class(d)

ans =

cell

>> d=[F{:}] %same can be done below

d =

1 1 1 0 3

1 1 1 0 4

>> d=[F{1},F{2},F{3}] %same as above

d =

1 1 1 0 3

1 1 1 0 4

>> e=[F{2:3}] %we can also extract subset of F

e =

0 3

0 4

>> e=[F{2},F{3}] %same as above

e =

0 3

0 4

>> celldisp(F)

F{1} =

1 1 1

F{2} =

F{3} =

>> [r,s,t]=deal(F{:}) %the function deal assigns the first variable r to first output F{1}, s to F{2}, t to F{3}

r =

1 1 1

s =

t =

>> [r,s,t]=deal(F{1},F{2},F{3}) %same as above

r =

1 1 1

s =

t =

>> [G{:}]=deal(a,b,c) %wrong way as we must specify dimensions of G as below

The left hand side has G{:} inside brackets, which requires that G be defined, so that the number of expected results

can be computed.

>> [G{1:3}]=deal(a,b,c)

G =

[2×3 double] [2×1 double] [2×1 double]

>> isequal(F,G) %they should be equal

ans =

2.5 Cell Functions

>> A(1,1)={[1:3;4:6;7:9]}, A(1,2)={2+3i}, A(2,1)={‘A character string’}, A(2,2)={12:-2:0}

A =

[3×3 double]

A =

[3×3 double] [2.0000 + 3.0000i]

A =

[3×3 double] [2.0000 + 3.0000i]

‘A character string’ []

A =

[3×3 double] [2.0000 + 3.0000i]

‘A character string’ [1×7 double]

>> cellfun(‘isreal’,A) %checks if a cell content is real or complex

ans =

1 0

1 1

>> cellfun(‘length’,A) %checks length of cells

ans =

3 1

18 7

>> cellfun(‘prodofsize’,A) %checks bumber of element in each cell

ans =

9 1

18 7

>> cellfun(‘isclass’,A,’char’) %check if a cell is a character string; there are other many cellfun functions

ans =

0 0

1 0

>> a=rand(3,6)

a =

0.8147 0.9134 0.2785 0.9649 0.9572 0.1419

0.9058 0.6324 0.5469 0.1576 0.4854 0.4218

0.1270 0.0975 0.9575 0.9706 0.8003 0.9157

>> c=num2cell(a) %builds a cell array c, that’s exactly as a note that a is not necessarily numbers as the name suggests

c =

[0.8147] [0.9134] [0.2785] [0.9649] [0.9572] [0.1419]

[0.9058] [0.6324] [0.5469] [0.1576] [0.4854] [0.4218]

[0.1270] [0.0975] [0.9575] [0.9706] [0.8003] [0.9157]

>> d=num2cell(a,1) %builds a cell array c from a by building a number of cells = number of columns of a and each column is a cell

d =

[3×1 double] [3×1 double] [3×1 double] [3×1 double] [3×1 double] [3×1 double]

>> e=num2cell(a,2) %builds a cell array c from a by building a number of cells = number of rows of a and each row is a cell

e =

[1×6 double]

2.6 Cell Arrays of Strings

>> s=’Hello World’

s =

Hello World

>> cs={‘Hi everyone’

‘Welcome to planet Earth’} %creates a cell array of strings

cs =

‘Hi everyone’

‘Hello World’

‘Welcome to planet Earth’

>> size(cs)

ans =

3 1

>> iscell(cs)

ans =

>> sa=char(cs) %builds a char string array from z cell array

sa =

Hi everyone

Hello World

Welcome to planet Earth

>> ischar(sa)

ans =

>> iscell(sa)

ans =

>> size(sa)

ans =

3 23

>> size(cs)

ans =

3 1

>> cst=cellstr(sa) %convert a string array into a cell array

cst =

‘Hi everyone’

‘Hello World’

‘Welcome to planet Earth’

>> iscell(cst)

ans =

>> size(cst)

ans =

3 1

>> isequal(cs,cst)

ans =

>> isequal(cs,sa)

ans =

2.7 Structure Creation

>> circle.radius=2.5, circle.center=[0 1], circle.linestyle=’–‘, circle.color=’red’ %we’re creating a structure here names circle that has 4 variables of different types; structures are pretty much like cell arrays, but they’re indexed using names not numbers

circle =

radius: 2.5000

circle =

radius: 2.5000

center: [0 1]

circle =

radius: 2.5000

center: [0 1]

linestyle: ‘–‘

circle =

radius: 2.5000

center: [0 1]

linestyle: ‘–‘

color: ‘red’

>> whos %see what MATLAB tells us about circle

Name Size Bytes Class Attributes

circle 1×1 530 struct

>> circle(2).radius=’sqrt(2)’, circle.center=[2.3 -12], circle.linestyle=’:’, circle.color=’green’ %we can create a second element in the structure circle, note that variable types don’t need to match; here, for instance, the radius variable is of type character

circle =

1×2 struct array with fields:

radius

center

linestyle

color

>> circle(2).radius=’sqrt(2)’, circle(2).center=[2.3 -12], circle(2).linestyle=’:’, circle(2).color=’green’ %we can create a second element in the structure circle, note that variable types don’t need to match; here, for instance, the radius variable is of type character

circle =

1×2 struct array with fields:

radius

center

linestyle

color

circle =

1×2 struct array with fields:

radius

center

linestyle

color

circle =

1×2 struct array with fields:

radius

center

linestyle

color

circle =

1×2 struct array with fields:

radius

center

linestyle

color

>> whos %now check size of circle

Name Size Bytes Class Attributes

circle 1×2 812 struct

>> circle(3).radius=25.4, circle(3).center=[-1 0], circle(3).linestyle=’-.’, circle(3).color=’blue’ %adding a 3rd element to circle structure, note that structure can take multidimensions, but it’s common to be of 1 dimension, compare adding 3 elements in a structure to how you would create this using arrays which is way more cumbersome

circle =

1×3 struct array with fields:

radius

center

linestyle

color

circle =

1×3 struct array with fields:

radius

center

linestyle

color

circle =

1×3 struct array with fields:

radius

center

linestyle

color

circle =

1×3 struct array with fields:

radius

center

linestyle

color

>> circle(1).filled=’yes’ %adding a 4th variable, that needs not be for all elements in the structure

circle =

1×3 struct array with fields:

radius

center

linestyle

color

filled

>> circle.filled %see how it looks like

ans =

yes

ans =

[]

ans =

[]

>> circle(2).filled=’no’, circle(3).filled=’yes’ %adding values to 2nd and 3rd elements

circle =

1×3 struct array with fields:

radius

center

linestyle

color

filled

circle =

1×3 struct array with fields:

radius

center

linestyle

color

filled

>> values1={2.5 ‘sqrt(2)’ 25.4}, values2={[0 1] [2.3 -1.2] [-1 0]}, values2={‘–‘ ‘:’ ‘-.’}, values4={‘red’ ‘green’ ‘blue’}, values5={‘yes’ ‘no’ ‘yes’} %creating cell arrays carrying values of variables

values1 =

[2.5000] ‘sqrt(2)’ [25.4000]

Values2 =

[1×2 double] [1×2 double] [1×2 double]

values2 =

‘–‘ ‘:’ ‘-.’

values4 =

‘red’ ‘green’ ‘blue’

values5 =

‘yes’ ‘no’ ‘yes’

>> values1={2.5 ‘sqrt(2)’ 25.4}, values2={[0 1] [2.3 -12] [-1 0]}, values3={‘–‘ ‘:’ ‘-.’}, values4={‘red’ ‘green’ ‘blue’}, values5={‘yes’ ‘no’ ‘yes’} %creating cell arrays carrying values of variables

values1 =

[2.5000] ‘sqrt(2)’ [25.4000]

values2 =

[1×2 double] [1×2 double] [1×2 double]

values3 =

‘–‘ ‘:’ ‘-.’

values4 =

‘red’ ‘green’ ‘blue’

values5 =

‘yes’ ‘no’ ‘yes’

>> CIRCLE=struct(‘radius’,values1,’center’,values2,’linestyle’,values3,’color’,values4,’filled’,values5) %we can create a structure using the function struct by providing cell arrays where values are and giving names to variables

CIRCLE =

1×3 struct array with fields:

radius

center

linestyle

color

filled

>> isequal(circle, CIRCLE) %check if they’re z same

ans =

2.8 Structure Manipulation

>> square.width=5, square.height=14, square.center=zeros(1,2), square.rotation=pi/4 %creating a new structure

square =

width: 5

square =

width: 5

height: 14

square =

width: 5

height: 14

center: [0 0]

square =

width: 5

height: 14

center: [0 0]

rotation: 0.7854

>> A=[circle CIRCLE] %concatenating circle and CIRCLE into a new structure; note that they must match in fields and names must be z same exactly

A =

1×6 struct array with fields:

radius

center

linestyle

color

filled

>> B=[circle square] %this shouldn’t work because of difference in fields

Error using horzcat

Number of fields in structure arrays being concatenated do not match. Concatenation of structure arrays requires that

these arrays have the same set of fields.

>> C=[circle(1:2) CIRCLE(3)] %we can also concatentaing subsets of structures

C =

1×3 struct array with fields:

radius

center

linestyle

color

filled

>> isequal(C, circle) %they share same fields and both have 3 elements

ans =

>> Aa=reshape(A,3,2) %reshape can be used here too

Aa =

3×2 struct array with fields:

radius

center

linestyle

color

filled

>> S=repmat(square,3,1) %repmat too

S =

3×1 struct array with fields:

width

height

center

rotation

2.9 Retrieving Structure Content

>> rad2=circle(2).radius

rad2 =

sqrt(2)

>> circle(1).radius

ans =

2.5000

>> area1=pi*circle(1).radius^2

area1 =

19.6350

>> circle(1).filled

ans =

yes

>> circle(1).filled(1)

ans =

>> circle(1).filled(2:end)

ans =

>> col=circle.color

col =

red

2.10 Comma-Separated Lists (Revisited)

>> a=ones(2,3),b=zeros(2,1), c=(3:4)’

a =

1 1 1

b =

c =

>> d=[a,b,c]

d =

1 1 1 0 3

1 1 1 0 4

>> d=cat(2,a,b,c)

d =

1 1 1 0 3

1 1 1 0 4

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

m =

n =

>> cent=cat(1,circle.center) %combines all fields of center in structure circle

cent =

0 1.0000

2.3000 -12.0000

-1.0000 0

>> cent=cat(1,circle(1).center,circle(2).center,circle(3).center)

cent =

0 1.0000

2.3000 -12.0000

-1.0000 0

>> some=cat(1,circle(2:end).center)

some =

2.3000 -12.0000

-1.0000 0

>> circle.color

ans =

red

ans =

green

ans =

blue

>> col=cat(1,circle.color) %they are not same lenght

Error using cat

Dimensions of matrices being concatenated are not consistent.

>> col=[circle.color]

col =

redgreenblue

>> col={circle.color}

col =

‘red’ ‘green’ ‘blue’

>> col=char(col)

col =

red

green

blue

>> [c1,c2,c3]=deal(circle.color) %assign the value of color in circle(1) to c1, circle(2) to c2 …

c1 =

red

c2 =

green

c3 =

blue

>> [r1,r2]=deal(circle([1 3]).radius) %assign the value of circle(1) radius to variable r1 and circle(3) to r2

r1 =

2.5000

r2 =

25.4000

>> [circle.radius]=deal(5,14,30) %assigns circle(1).radius to 5, circle(2).radius to 14 …

circle =

1×3 struct array with fields:

radius

center

linestyle

color

filled

>> circle.radius

ans =

>> [triangle(:).type]=deal(‘right’,’isosceles’,’unknown’) %size of triangle is not given, therefore an error is thrown

Error using deal (line 38)

The number of outputs should match the number of inputs.

>> [triangle(1:3).type]=deal(‘right’,’isosceles’,’unknown’) %size is given, so it works

triangle =

1×3 struct array with fields:

type

>> triangle.type

ans =

right

ans =

isosceles

ans =

unknown

2.11 Structure Functions

>> circle

circle =

1×3 struct array with fields:

radius

center

linestyle

color

filled

>> square

square =

width: 5

height: 14

center: [0 0]

rotation: 0.7854

>> fieldnames(circle) %lists field names of circle

ans =

‘radius’

‘center’

‘linestyle’

‘color’

‘filled’

>> isfield(circle,’color’) %check if color is a field name

ans =

>> isfield(circle,’width’)

ans =

>> class(square)

ans =

struct

>> isstruct(circle) %check if circle is a structure

ans =

>> isstruct(r1)

ans =

>> fnames=fieldnames(circle) %assigns names of fields of circle structure to the array fnames

fnames =

‘radius’

‘center’

‘linestyle’

‘color’

‘filled’

>> circle2=rmfield(circle,fnames{5}) %rmfields removes fields, here it removes the 5th field (filled)

circle2 =

1×3 struct array with fields:

radius

center

linestyle

color

>> r1=getfield(circle,{1},fnames{1}) %getfield brings the value in a field, remember structures are indexed with names

r1 =

>> r2=getfield(circle,{3},fnames{1})

r2 =

>> circle3=setfield(circle,{3},fnames{1},r2) %setfield set the value of a field to a certain value

circle3 =

1×3 struct array with fields:

radius

center

linestyle

color

filled

>> circle3(3).radius

ans =

4. Integration and Differentiation

4.1 Integration

>> x=-1:.17:2;

>> y=humps(x);

>> area=trapz(x,y) %finds approximate value for integration by estimating area under curve using trapezodial approximation

area =

25.9174

>> x=linspace(-1,2,100);

>> y=humps(x);

>> fomrat long

Undefined function ‘fomrat’ for input arguments of type ‘char’.

Did you mean:

>> format long

>> area=trapz(x,y) %here more accurate result is obtained because we have a finer discretization

area =

26.344731195245952

>> x=linspace(-1,2,100);

>> y=humps(x);

>> z=cumtrapz(x,y); %cumtrapz function finds cumulative integrals

>> size(z)

ans =

1 100

>> plotyy(x,y,x,z)

>> grid on

>> z(end) %finds value of integration

ans =

26.344731195245959

>> quad(@humps,-1,2) %quad function optimizes the width of trapezoids to best fit the function characteristics

ans =

26.344960501201232

>> quadl(@humps,-1,2) %quad and quadl give more accurate results; quadl is more rigorous

ans =

26.344960471378968

********************************M-FILE*************************

function z=myfun(x,y)

z=sin(x).*cos(y)+1;

********************************M-FILE*************************

>> x=linspace(0,pi,20);

>> y=linspace(-pi,pi,20);

>> [xx,yy]=meshgrid(x,y);

>> zz=myfun(xx,yy);

>> mesh(xx,yy,zz)

>> area=dblquad(‘myfun’,0,pi,-pi,pi) %finding double integration

area =

19.739208806091021

>> relerror=(area-2*pi^2)/(2*pi^2) %finding relative error

relerror =

1.981996941074027e-10

4.2 Differentiation

>> x=0:0.1:1 ;

>> y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2] ;

>> p=polyfit(x,y,2) %find a polynomial of degree 2 to fit the data

p =

-9.8108 20.1293 -0.0317

>> xi=linspace(0,1,100);

>> yi=polyval(p,xi); %evaluate the polynomial over the xi interval

>> plot(x,y,’-o’,xi,yi,’–‘)

>> pd=polyder(p) %find polynomial derivative

pd =

-19.6217 20.1293

>> dyp=polyval(pd,x); %evaluate the derivative polynomial

>> dy=diff(y)./diff(x); % evaluate derivative using forward difference and divsion (recall derivative = dy/dx)

>> xd=x(1:end-1); %generate new array for dy as number of dy elements is less than those of y by 1 (difference)

>> plot(xd,dy,x,dyp,’:’)

>> x=linspace(0,2*pi);

>> y=sin(x);

>> dy=diff(y)/(x(2)-x(1)); %here x is linearly spaced so no need for diff(x)

>> xd=x(2:end); %here we are using backward difference

>> plot(x,y,xd,dy)

>> max(abs(cos(xd)-dy)) %finding maximum error

ans =

0.0317

>> dy=(y(3:end)-y(1:end-2))/(x(3)-x(1)); %now, central differencing

>> xd=x(2:end-1);

>> max(abs(cos(xd)-dy)) %finding maximum error in central difference case

ans =

6.7086e-04

>> [x,y,z]=peaks(20); %tabualted data as z=f(x,y)

>> dx=x(1,2)-x(1,1) ;

>> dy=y(2,1)-y(1,1) ;

>> [dzdx,dzdy]=gradient(z,dx,dy); %find dz/dx and dz/dy (partial derivatives)

>> contour(x,y,z), hold on, quiver(x,y,dzdx,dzdy) %quiver plots slope at each point with length scaled to value and arrow head to direction

>> [x,y,z]=peaks; %tabualted data as z=f(x,y)

>> dx=x(1,2)-x(1,1);

>> dy=y(2,1)-y(1,1) ;

>> L=del2(z,dx,dy); %del2 calculates the Laplacian for z as a measure for surface curvature

>> surf(x,y,z,abs(L)), shading interp %colors represent the degree of curvature

5. Differential Equations

5.1 IVP Format

The standard format is:

Our sample equation is:

Set y₁ = x and y₂ = dx/dt

So, the 2 equations are:

5.2 ODE Suite Solvers

>> help ode23

ode23 Solve non-stiff differential equations, low order method.

[TOUT,YOUT] = ode23(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates

the system of differential equations y’ = f(t,y) from time T0 to TFINAL

with initial conditions Y0. ODEFUN is a function handle. For a scalar T

and a vector Y, ODEFUN(T,Y) must return a column vector corresponding

to f(t,y). Each row in the solution array YOUT corresponds to a time

returned in the column vector TOUT. To obtain solutions at specific

times T0,T1,…,TFINAL (all increasing or all decreasing), use TSPAN =

[T0 T1 … TFINAL].

[TOUT,YOUT] = ode23(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default

integration properties replaced by values in OPTIONS, an argument created

with the ODESET function. See ODESET for details. Commonly used options

are scalar relative error tolerance ‘RelTol’ (1e-3 by default) and vector

of absolute error tolerances ‘AbsTol’ (all components 1e-6 by default).

If certain components of the solution must be non-negative, use

ODESET to set the ‘NonNegative’ property to the indices of these

components.

ode23 can solve problems M(t,y)*y’ = f(t,y) with mass matrix M that is

nonsingular. Use ODESET to set the ‘Mass’ property to a function handle

MASS if MASS(T,Y) returns the value of the mass matrix. If the mass matrix

is constant, the matrix can be used as the value of the ‘Mass’ option. If

the mass matrix does not depend on the state variable Y and the function

MASS is to be called with one input argument T, set ‘MStateDependence’ to

‘none’. ODE15S and ODE23T can solve problems with singular mass matrices.

[TOUT,YOUT,TE,YE,IE] = ode23(ODEFUN,TSPAN,Y0,OPTIONS) with the ‘Events’

property in OPTIONS set to a function handle EVENTS, solves as above

while also finding where functions of (T,Y), called event functions,

are zero. For each function you specify whether the integration is

to terminate at a zero and whether the direction of the zero crossing

matters. These are the three column vectors returned by EVENTS:

[VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function:

VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration

is to terminate at a zero of this event function and 0 otherwise.

DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only

zeros where the event function is increasing, and -1 if only zeros where

the event function is decreasing. Output TE is a column vector of times

at which events occur. Rows of YE are the corresponding solutions, and

indices in vector IE specify which event occurred.

SOL = ode23(ODEFUN,[T0 TFINAL],Y0…) returns a structure that can be

used with DEVAL to evaluate the solution or its first derivative at

any point between T0 and TFINAL. The steps chosen by ode23 are returned

in a row vector SOL.x. For each I, the column SOL.y(:,I) contains

the solution at SOL.x(I). If events were detected, SOL.xe is a row vector

of points at which events occurred. Columns of SOL.ye are the corresponding

solutions, and indices in vector SOL.ie specify which event occurred.

Example

[t,y]=ode23(@vdp1,[0 20],[2 0]);

plot(t,y(:,1));

solves the system y’ = vdp1(t,y), using the default relative error

tolerance 1e-3 and the default absolute tolerance of 1e-6 for each

component, and plots the first component of the solution.

Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):

float: double, single

See also ode45, ode113, ode15s, ode23s, ode23t, ode23tb, ode15i,

odeset, odeplot, odephas2, odephas3, odeprint, deval,

odeexamples, rigidode, ballode, orbitode, function_handle.

5.3 Basic Use

********************************M-FILE*************************

function ydot=vdpol(t,y,mu)

if nargin<3

mu=2;

end

ydot=[y(2); mu*(1-y(1)^2)*y(2)-y(1)];

********************************M-FILE*************************

>> tspan=[0 20]; %time span

>> yo=[2; 0]; %initial condiitons (note it must be a column)

>> [t,y]=ode45(@vdpol,tspan,yo); %solving the set of 2 differential equations using the ode45 solver; output will be stored in t and y arrays; function name (vdpol) must be given that contains the differential equations and also the time span and initial conditions must be given

>> plot(t,y(:,1),t,y(:,2),’–‘) %plot t-vs-y1 and t-vs-y2

>> mu=10

mu =

>> ode45(@vdpol,tspan,yo,mu); %solving using ode45 and with mu=10, note that no output variables are set, so it just plots the solution

>> ode45(@vdpol,tspan,yo,[],mu); %solving using ode45 and with mu=10, note that no output variables are set, so it just plots the solution

5.4 Setting Options

Options

>> help odeset

odeset Create/alter ODE OPTIONS structure.

OPTIONS = odeset(‘NAME1′,VALUE1,’NAME2′,VALUE2,…) creates an integrator

options structure OPTIONS in which the named properties have the

specified values. Any unspecified properties have default values. It is

sufficient to type only the leading characters that uniquely identify the

property. Case is ignored for property names.

OPTIONS = odeset(OLDOPTS,’NAME1’,VALUE1,…) alters an existing options

structure OLDOPTS.

OPTIONS = odeset(OLDOPTS,NEWOPTS) combines an existing options structure

OLDOPTS with a new options structure NEWOPTS. Any new properties

overwrite corresponding old properties.

odeset with no input arguments displays all property names and their

possible values.

odeset PROPERTIES

RelTol – Relative error tolerance [ positive scalar {1e-3} ]

This scalar applies to all components of the solution vector, and

defaults to 1e-3 (0.1% accuracy) in all solvers. The estimated error in

each integration step satisfies e(i) <= max(RelTol*abs(y(i)),AbsTol(i)).

AbsTol – Absolute error tolerance [ positive scalar or vector {1e-6} ]

A scalar tolerance applies to all components of the solution vector.

Elements of a vector of tolerances apply to corresponding components of

the solution vector. AbsTol defaults to 1e-6 in all solvers. See RelTol.

NormControl – Control error relative to norm of solution [ on | {off} ]

Set this property ‘on’ to request that the solvers control the error in

each integration step with norm(e) <= max(RelTol*norm(y),AbsTol). By

default the solvers use a more stringent component-wise error control.

Refine – Output refinement factor [ positive integer ]

This property increases the number of output points by the specified

factor producing smoother output. Refine defaults to 1 in all solvers

except ODE45, where it is 4. Refine does not apply if length(TSPAN) > 2

or the ODE solver returns the solution as a structure.

OutputFcn – Installable output function [ function_handle ]

This output function is called by the solver after each time step. When

a solver is called with no output arguments, OutputFcn defaults to

@odeplot. Otherwise, OutputFcn defaults to [].

OutputSel – Output selection indices [ vector of integers ]

This vector of indices specifies which components of the solution vector

are passed to the OutputFcn. OutputSel defaults to all components.

Stats – Display computational cost statistics [ on | {off} ]

Jacobian – Jacobian function [ function_handle | constant matrix ]

Set this property to @FJac if FJac(t,y) returns dF/dy, or to

the constant value of dF/dy.

For ODE15I solving F(t,y,y’) = 0, set this property to @FJac if

[dFdy, dFdyp] = FJac(t,y,yp), or to a cell array of constant

values {dF/dy,dF/dyp}.

JPattern – Jacobian sparsity pattern [ sparse matrix ]

Set this property to a sparse matrix S with S(i,j) = 1 if component i of

F(t,y) depends on component j of y, and 0 otherwise.

For ODE15I solving F(t,y,y’) = 0, set this property to

{dFdyPattern,dFdypPattern}, the sparsity patterns of dF/dy and

dF/dy’, respectively.

Vectorized – Vectorized ODE function [ on | {off} ]

Set this property ‘on’ if the ODE function F is coded so that

F(t,[y1 y2 …]) returns [F(t,y1) F(t,y2) …].

For ODE15I solving F(t,y,y’) = 0, set this property to

{yVect,ypVect}. Setting yVect ‘on’ indicates that

F(t,[y1 y2 …],yp) returns [F(t,y1,yp) F(t,y2,yp) …].

Setting ypVect ‘on’ indicates that F(t,y,[yp1 yp2 …])

returns [F(t,y,yp1) F(t,y,yp2) …].

Events – Locate events [ function_handle ]

To detect events, set this property to the event function.

Mass – Mass matrix [ constant matrix | function_handle ]

For problems M*y’ = f(t,y) set this property to the value of the constant

mass matrix. For problems with time- or state-dependent mass matrices,

set this property to a function that evaluates the mass matrix.

MStateDependence – Dependence of the mass matrix on y [ none | {weak} | strong ]

Set this property to ‘none’ for problems M(t)*y’ = F(t,y). Both ‘weak’ and

‘strong’ indicate M(t,y), but ‘weak’ will result in implicit solvers

using approximations when solving algebraic equations.

MassSingular – Mass matrix is singular [ yes | no | {maybe} ]

Set this property to ‘no’ if the mass matrix is not singular.

MvPattern – dMv/dy sparsity pattern [ sparse matrix ]

Set this property to a sparse matrix S with S(i,j) = 1 if for any k, the

(i,k) component of M(t,y) depends on component j of y, and 0 otherwise.

InitialSlope – Consistent initial slope yp0 [ vector ]

yp0 satisfies M(t0,y0)*yp0 = F(t0,y0).

InitialStep – Suggested initial step size [ positive scalar ]

The solver will try this first. By default the solvers determine an

initial step size automatically.

MaxStep – Upper bound on step size [ positive scalar ]

MaxStep defaults to one-tenth of the tspan interval in all solvers.

NonNegative – Non-negative solution components [ vector of integers ]

This vector of indices specifies which components of the

solution vector must be non-negative. NonNegative defaults to [].

This property is not available in ODE23S, ODE15I. In ODE15S,

ODE23T, and ODE23TB, the property is not available for problems

where there is a mass matrix.

BDF – Use Backward Differentiation Formulas in ODE15S [ on | {off} ]

This property specifies whether the Backward Differentiation Formulas

(Gear’s methods) are to be used in ODE15S instead of the default

Numerical Differentiation Formulas.

MaxOrder – Maximum order of ODE15S and ODE15I [ 1 | 2 | 3 | 4 | {5} ]

See also odeget, ode45, ode23, ode113, ode15i, ode15s, ode23s, ode23t, ode23tb,

function_handle.

>> odeset

AbsTol: [ positive scalar or vector {1e-6} ]

RelTol: [ positive scalar {1e-3} ]

NormControl: [ on | {off} ]

NonNegative: [ vector of integers ]

OutputFcn: [ function_handle ]

OutputSel: [ vector of integers ]

Refine: [ positive integer ]

Stats: [ on | {off} ]

InitialStep: [ positive scalar ]

MaxStep: [ positive scalar ]

BDF: [ on | {off} ]

MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]

Jacobian: [ matrix | function_handle ]

JPattern: [ sparse matrix ]

Vectorized: [ on | {off} ]

Mass: [ matrix | function_handle ]

MStateDependence: [ none | {weak} | strong ]

MvPattern: [ sparse matrix ]

MassSingular: [ yes | no | {maybe} ]

InitialSlope: [ vector ]

Events: [ function_handle ]

Setting tolerance

>> tspan=[0 20]; %time span

>> yo=[2; 0]; %initial condiitons

>> mu=10;

>> options=odeset(‘AbsTol’,1e-12,’RelTol’,1e-6)

options =

AbsTol: 1.0000e-12

BDF: []

Events: []

InitialStep: []

Jacobian: []

JConstant: []

JPattern: []

Mass: []

MassSingular: []

MaxOrder: []

MaxStep: []

NonNegative: []

NormControl: []

OutputFcn: []

OutputSel: []

Refine: []

RelTol: 1.0000e-06

Stats: []

Vectorized: []

MStateDependence: []

MvPattern: []

InitialSlope: []

>> [t,y]=ode45(@vdpol,tspan,yo,[],mu); %default tolerances; the [] tells MATLAB to use default options

>> length(t)

ans =

593

>> [t,y]=ode45(@vdpol,tspan,yo,options,mu); %tight tolerances (set in options)

>> length(t)

ans =

1689

>> [t,y]=ode15s(@vdpol,tspan,yo,[],mu); %different the solver to the stiff ode15s

>> length(t)

ans =

232

>> [t,y]=ode15s(@vdpol,tspan,yo,options,mu); %ode15s solver with tight tolerances

>> length(t)

ans =

651

Jacobian for solvers

********************************M-FILE*************************

function jac=vdpoljac(t,y,mu)

if nargin<3

mu=2;

end

jac= [0 1

(-2*mu*y(1)*y(2)-1) (mu*(1-y(1)^2))];

********************************M-FILE*************************

>> options=odeset(options,’Jacobian’,@vdpoljac) %setting the Jacobian to the function vdpoljac; giving MATLAB an analytical Jacobian in cases of stiff solvers is always recommended

options =

AbsTol: 1.0000e-12

BDF: []

Events: []

InitialStep: []

Jacobian: @vdpoljac

JConstant: []

JPattern: []

Mass: []

MassSingular: []

MaxOrder: []

MaxStep: []

NonNegative: []

NormControl: []

OutputFcn: []

OutputSel: []

Refine: []

RelTol: 1.0000e-06

Stats: []

Vectorized: []

MStateDependence: []

MvPattern: []

InitialSlope: []

>> [t,y]=ode15s(@vdpol,tspan,yo,options,mu); %ode15s solver with Jacobian given

>> length(t) %even though it is more than previous, execution speed increases

ans =

670

Vectorized ODE and Refine property

********************************M-FILE*************************

function ydot=vdpolv(t,y,mu)

if nargin<3

mu=2;

end

ydot=[y(2,:); mu*(1-y(1,:)^2)*y(2,:)-y(1,:)];

********************************M-FILE*************************

>> options=odeset(‘Refine’,1) %refine sets how many intermediate points to interpolate the solution at within each integration step

options =

AbsTol: []

BDF: []

Events: []

InitialStep: []

Jacobian: []

JConstant: []

JPattern: []

Mass: []

MassSingular: []

MaxOrder: []

MaxStep: []

NonNegative: []

NormControl: []

OutputFcn: []

OutputSel: []

Refine: 1

RelTol: []

Stats: []

Vectorized: []

MStateDependence: []

MvPattern: []

InitialSlope: []

>> [t,y]=ode45(@vdpolv,tspan,yo,options,mu); %solving with vectorized ODE file

>> length(t)

ans =

149

Events property

********************************M-FILE*************************

function [value,isterminal,direction]=vdpolevents(t,y,mu)

value(1)=abs(y(2))-2;

isterminal(1)=0;

direction(1)=0;

********************************M-FILE*************************

>> mu=2;

>> options=odeset(‘Events’,@vdpolevents); %calls the vdpolevents function that tells us when abs(y(2)) = 2

>> [t,y,te,ye]=ode45(@vdpolv,tspan,yo,options,mu); %solving with events parameter; note the added te and ye for events

>> plot(t,y,te,ye(:,2),’o’)

5.5 BVPs and PDEs

>> help bvp4c

bvp4c Solve boundary value problems for ODEs by collocation.

SOL = bvp4c(ODEFUN,BCFUN,SOLINIT) integrates a system of ordinary

differential equations of the form y’ = f(x,y) on the interval [a,b],

subject to general two-point boundary conditions of the form

bc(y(a),y(b)) = 0. ODEFUN and BCFUN are function handles. For a scalar X

and a column vector Y, ODEFUN(X,Y) must return a column vector representing

f(x,y). For column vectors YA and YB, BCFUN(YA,YB) must return a column

vector representing bc(y(a),y(b)). SOLINIT is a structure with fields

x — ordered nodes of the initial mesh with

SOLINIT.x(1) = a, SOLINIT.x(end) = b

y — initial guess for the solution with SOLINIT.y(:,i)

a guess for y(x(i)), the solution at the node SOLINIT.x(i)

bvp4c produces a solution that is continuous on [a,b] and has a

continuous first derivative there. The solution is evaluated at points

XINT using the output SOL of bvp4c and the function DEVAL:

YINT = DEVAL(SOL,XINT). The output SOL is a structure with

SOL.solver — ‘bvp4c’

SOL.x — mesh selected by bvp4c

SOL.y — approximation to y(x) at the mesh points of SOL.x

SOL.yp — approximation to y'(x) at the mesh points of SOL.x

SOL.stats — computational cost statistics (also displayed when

the ‘Stats’ option is set with BVPSET).

SOL = bvp4c(ODEFUN,BCFUN,SOLINIT,OPTIONS) solves as above with default

parameters replaced by values in OPTIONS, a structure created with the

BVPSET function. To reduce the run time greatly, use OPTIONS to supply

a function for evaluating the Jacobian and/or vectorize ODEFUN.

See BVPSET for details and SHOCKBVP for an example that does both.

Some boundary value problems involve a vector of unknown parameters p

that must be computed along with y(x):

y’ = f(x,y,p)

0 = bc(y(a),y(b),p)

For such problems the field SOLINIT.parameters is used to provide a guess

for the unknown parameters. On output the parameters found are returned

in the field SOL.parameters. The solution SOL of a problem with one set

of parameter values can be used as SOLINIT for another set. Difficult BVPs

may be solved by continuation: start with parameter values for which you can

get a solution, and use it as a guess for the solution of a problem with

parameters closer to the ones you want. Repeat until you solve the BVP

for the parameters you want.

The function BVPINIT forms the guess structure in the most common

situations: SOLINIT = BVPINIT(X,YINIT) forms the guess for an initial

mesh X as described for SOLINIT.x, and YINIT either a constant vector

guess for the solution or a function handle. If YINIT is a function handle

then for a scalar X, YINIT(X) must return a column vector, a guess for

the solution at point x in [a,b]. If the problem involves unknown parameters

SOLINIT = BVPINIT(X,YINIT,PARAMS) forms the guess with the vector PARAMS of

guesses for the unknown parameters.

bvp4c solves a class of singular BVPs, including problems with

unknown parameters p, of the form

y’ = S*y/x + f(x,y,p)

0 = bc(y(0),y(b),p)

The interval is required to be [0, b] with b > 0.

Often such problems arise when computing a smooth solution of

ODEs that result from PDEs because of cylindrical or spherical

symmetry. For singular problems the (constant) matrix S is

specified as the value of the ‘SingularTerm’ option of BVPSET,

and ODEFUN evaluates only f(x,y,p). The boundary conditions

must be consistent with the necessary condition S*y(0) = 0 and

the initial guess should satisfy this condition.

bvp4c can solve multipoint boundary value problems. For such problems

there are boundary conditions at points in [a,b]. Generally these points

represent interfaces and provide a natural division of [a,b] into regions.

bvp4c enumerates the regions from left to right (from a to b), with indices

starting from 1. In region k, bvp4c evaluates the derivative as

YP = ODEFUN(X,Y,K). In the boundary conditions function,

BCFUN(YLEFT,YRIGHT), YLEFT(:,K) is the solution at the ‘left’ boundary

of region k and similarly for YRIGHT(:,K). When an initial guess is

created with BVPINIT(XINIT,YINIT), XINIT must have double entries for

each interface point. If YINIT is a function handle, BVPINIT calls

Y = YINIT(X,K) to get an initial guess for the solution at X in region k.

In the solution structure SOL returned by bvp4c, SOL.x has double entries

for each interface point. The corresponding columns of SOL.y contain

the ‘left’ and ‘right’ solution at the interface, respectively.

See THREEBVP for an example of solving a three-point BVP.

Example

solinit = bvpinit([0 1 2 3 4],[1 0]);

sol = bvp4c(@twoode,@twobc,solinit);

solve a BVP on the interval [0,4] with differential equations and

boundary conditions computed by functions twoode and twobc, respectively.

This example uses [0 1 2 3 4] as an initial mesh, and [1 0] as an initial

approximation of the solution components at the mesh points.

xint = linspace(0,4);

yint = deval(sol,xint);

evaluate the solution at 100 equally spaced points in [0 4]. The first

component of the solution is then plotted with

plot(xint,yint(1,:));

For more examples see TWOBVP, FSBVP, SHOCKBVP, MAT4BVP, EMDENBVP, THREEBVP.

See also bvp5c, bvpset, bvpget, bvpinit, bvpxtend, deval, function_handle.

>> help pdepe

pdepe Solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D.

SOL = pdepe(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN) solves initial-boundary

value problems for small systems of parabolic and elliptic PDEs in one

space variable x and time t to modest accuracy. There are npde unknown

solution components that satisfy a system of npde equations of the form

c(x,t,u,Du/Dx) * Du/Dt = x^(-m) * D(x^m * f(x,t,u,Du/Dx))/Dx + s(x,t,u,Du/Dx)

Here f(x,t,u,Du/Dx) is a flux and s(x,t,u,Du/Dx) is a source term. m must

be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry,

respectively. The coupling of the partial derivatives with respect to

time is restricted to multiplication by a diagonal matrix c(x,t,u,Du/Dx).

The diagonal elements of c are either identically zero or positive.

An entry that is identically zero corresponds to an elliptic equation and

otherwise to a parabolic equation. There must be at least one parabolic

equation. An entry of c corresponding to a parabolic equation is permitted

to vanish at isolated values of x provided they are included in the mesh

XMESH, and in particular, is always allowed to vanish at the ends of the

interval. The PDEs hold for t0 <= t <= tf and a <= x <= b. The interval

[a,b] must be finite. If m > 0, it is required that 0 <= a. The solution

components are to have known values at the initial time t = t0, the

initial conditions. The solution components are to satisfy boundary

conditions at x=a and x=b for all t of the form

p(x,t,u) + q(x,t) * f(x,t,u,Du/Dx) = 0

q(x,t) is a diagonal matrix. The diagonal elements of q must be either

identically zero or never zero. Note that the boundary conditions are

expressed in terms of the flux rather than Du/Dx. Also, of the two

coefficients, only p can depend on u.

The input argument M defines the symmetry of the problem. PDEFUN, ICFUN,

and BCFUN are function handles.

[C,F,S] = PDEFUN(X,T,U,DUDX) evaluates the quantities defining the

differential equation. The input arguments are scalars X and T and

vectors U and DUDX that approximate the solution and its partial

derivative with respect to x, respectively. PDEFUN returns column

vectors: C (containing the diagonal of the matrix c(x,t,u,Dx/Du)),

F, and S (representing the flux and source term, respectively).

U = ICFUN(X) evaluates the initial conditions. For a scalar X, ICFUN

must return a column vector, corresponding to the initial values of

the solution components at X.

[PL,QL,PR,QR] = BCFUN(XL,UL,XR,UR,T) evaluates the components of the

boundary conditions at time T. XL and XR are scalars representing the

left and right boundary points. UL and UR are column vectors with the

solution at the left and right boundary, respectively. PL and QL are

column vectors corresponding to p and the diagonal of q, evaluated at

the left boundary, similarly PR and QR correspond to the right boundary.

When m > 0 and a = 0, boundedness of the solution near x = 0 requires

that the flux f vanish at a = 0. pdepe imposes this boundary condition

automatically.

pdepe returns values of the solution on a mesh provided as the input

array XMESH. The entries of XMESH must satisfy

a = XMESH(1) < XMESH(2) < … < XMESH(NX) = b

for some NX >= 3. Discontinuities in c and/or s due to material

interfaces are permitted if the problem requires the flux f to be

continuous at the interfaces and a mesh point is placed at each

interface. The ODEs resulting from discretization in space are integrated

to obtain approximate solutions at times specified in the input array

TSPAN. The entries of TSPAN must satisfy

t0 = TSPAN(1) < TSPAN(2) < … < TSPAN(NT) = tf

for some NT >= 3. The arrays XMESH and TSPAN do not play the same roles

in pdepe: The time integration is done with an ODE solver that selects

both the time step and formula dynamically. The cost depends weakly on

the length of TSPAN. Second order approximations to the solution are made

on the mesh specified in XMESH. Generally it is best to use closely

spaced points where the solution changes rapidly. pdepe does not select

the mesh in x automatically like it does in t; you must choose an

appropriate fixed mesh yourself. The discretization takes into account

the coordinate singularity at x = 0 when m > 0, so it is not necessary to

use a fine mesh near x = 0 for this reason. The cost depends strongly on

the length of XMESH.

The solution is returned as a multidimensional array SOL. UI = SOL(:,:,i)

is an approximation to component i of the solution vector u for

i = 1:npde. The entry UI(j,k) = SOL(j,k,i) approximates UI at

(t,x) = (TSPAN(j),XMESH(k)).

SOL = pdepe(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN,OPTIONS) solves as above

with default integration parameters replaced by values in OPTIONS, an

argument created with the ODESET function. Only some of the options of

the underlying ODE solver are available in pdepe – RelTol, AbsTol,

NormControl, InitialStep, and MaxStep. See ODESET for details.

[SOL,TSOL,SOLE,TE,IE] = pdepe(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN,OPTIONS)

with the ‘Events’ property in OPTIONS set to a function handle EVENTS,

solves as above while also finding where event functions g(t,u(x,t))

are zero. For each function you specify whether the integration is to

terminate at a zero and whether the direction of the zero crossing

matters. These are the three column vectors returned by EVENTS:

[VALUE,ISTERMINAL,DIRECTION] = EVENTS(M,T,XMESH,UMESH).

XMESH contains the spatial mesh and UMESH is the solution at the mesh

points. Use PDEVAL to evaluate the solution between mesh points.

For the I-th event function: VALUE(I) is the value of the function,

ISTERMINAL(I) = 1 if the integration is to terminate at a zero of this

event function and 0 otherwise. DIRECTION(I) = 0 if all zeros are to be

computed (the default), +1 if only zeros where the event function is

increasing, and -1 if only zeros where the event function is decreasing.

Output TSOL is a column vector of times specified in TSPAN, prior to

first terminal event. SOL(j,:,:) is the solution at T(j). TE is a vector

of times at which events occur. SOLE(j,:,:) is the solution at TE(j) and

indices in vector IE specify which event occurred.

If UI = SOL(j,:,i) approximates component i of the solution at time TSPAN(j)

and mesh points XMESH, PDEVAL evaluates the approximation and its partial

derivative Dui/Dx at the array of points XOUT and returns them in UOUT

and DUOUTDX: [UOUT,DUOUTDX] = PDEVAL(M,XMESH,UI,XOUT)

NOTE that the partial derivative Dui/Dx is evaluated here rather than the

flux. The flux is continuous, but at a material interface the partial

derivative may have a jump.

Example

x = linspace(0,1,20);

t = [0 0.5 1 1.5 2];

sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

solve the problem defined by the function pdex1pde with the initial and

boundary conditions provided by the functions pdex1ic and pdex1bc,

respectively. The solution is obtained on a spatial mesh of 20 equally

spaced points in [0 1] and it is returned at times t = [0 0.5 1 1.5 2].

Often a good way to study a solution is plot it as a surface and use

Rotate 3D. The first unknown, u1, is extracted from sol and plotted with

u1 = sol(:,:,1);

surf(x,t,u1);

PDEX1 shows how this problem can be coded using subfunctions. For more

examples see PDEX2, PDEX3, PDEX4, PDEX5. The examples can be read

separately, but read in order they form a mini-tutorial on using pdepe.

See also pdeval, ode15s, odeset, odeget, function_handle.

Reference page in Help browser

doc pdepe

6. Fourier Analysis

6.1 Discrete Fourier Transform

>> help fft

fft Discrete Fourier transform.

fft(X) is the discrete Fourier transform (DFT) of vector X. For

matrices, the fft operation is applied to each column. For N-D

arrays, the fft operation operates on the first non-singleton

dimension.

fft(X,N) is the N-point fft, padded with zeros if X has less

than N points and truncated if it has more.

fft(X,[],DIM) or fft(X,N,DIM) applies the fft operation across the

dimension DIM.

For length N input vector x, the DFT is a length N vector X,

with elements

X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.

n=1

The inverse DFT (computed by IFFT) is given by

x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.

k=1

See also fft2, fftn, fftshift, fftw, ifft, ifft2, ifftn.

Overloaded methods:

uint8/fft

uint16/fft

gf/fft

codistributed/fft

gpuArray/fft

qfft/fft

iddata/fft

>> N=128; %N is the length of the interval we computer over; always choose a number to be a power of 2 for fast computing

>> t=linspace(0,3,N); %time interval we compute over

>> f=2*exp(-3*t); %our original function in time domain

>> Ts=t(2)-t(1) %sampling period

Ts =

0.0236

>> Ws=2*pi/Ts %sampling frequency in rad/sec

Ws =

265.9882

>> W=Ws*(-N/2:(N/2)-1)/N; %frequency axis; note that W here is an array

>> F=fft(f); %compute Fourier transform using fft function

>> Fc=fftshift(F)*Ts; %now that’s shifting and scaling

>> Fa=2./(3+j*W); %that’s Fourier transform analytically

>> plot(W,abs(Fa),W,abs(Fc),’o’)

>> help conv

conv Convolution and polynomial multiplication.

C = conv(A, B) convolves vectors A and B. The resulting vector is

length MAX([LENGTH(A)+LENGTH(B)-1,LENGTH(A),LENGTH(B)]). If A and B are

vectors of polynomial coefficients, convolving them is equivalent to

multiplying the two polynomials.

C = conv(A, B, SHAPE) returns a subsection of the convolution with size

specified by SHAPE:

‘full’ – (default) returns the full convolution,

‘same’ – returns the central part of the convolution

that is the same size as A.

‘valid’ – returns only those parts of the convolution

that are computed without the zero-padded edges.

LENGTH(C)is MAX(LENGTH(A)-MAX(0,LENGTH(B)-1),0).

Class support for inputs A,B:

float: double, single

See also deconv, conv2, convn, filter and,

in the signal Processing Toolbox, xcorr, convmtx.

6.2 Fourier Series

********************************M-FILE*************************

function f=sawtooth(t,To)

f=10*rem(t,To)/To;

f(f==0 | f==10)=5;

********************************M-FILE*************************

>> N=25; %number of harmonics

>> To=0.2; %period

>> n=2*N; %number of terms to compute in the discrete Fourier transform is double that of harmonics, since it computes both positive and negative harmonics

>> t=linspace(0,To,n+1); %the function should be evaluated at n points over 1 period such that (n+1)th point be 1 period away from 1st point

>> t(end)=[]; %delete last opint (undesired)

>> f=sawtooth(t,To); %evaluate sawtooth function

>> Fn=fft(f); %Fourier transform using fft function

>> Fn=[conj(Fn(N+1)) Fn(N+2:end) Fn(1:N+1)]; %re-arranging the terms so that Fn(1) is F(-25) …

>> Fn=Fn/n; %scaling

>> A0=Fn(N+1) %note that Fn is an array of complex exponential Fourier series coefficients; Ao is trignometric representation first coefficients

A0 =

>> An=2*real(Fn(N+2:end)) %cosine terms

An =

1.0e-15 *

Columns 1 through 12

-0.4086 0.0622 -0.1870 0.2130 0.0439 0.1454 0.0188 -0.0664 0.0378 0 -0.1292 -0.0664

Columns 13 through 24

-0.0578 -0.0402 -0.1150 0.0015 0.0485 0.0131 -0.2404 0 0.0299 0.0131 0.2220 0.3908

Column 25

0.1421

>> Bn=-2*imag(Fn(N+2:end)) %sine terms

Bn =

Columns 1 through 11

-3.1789 -1.5832 -1.0484 -0.7789 -0.6155 -0.5051 -0.4250 -0.3638 -0.3151 -0.2753 -0.2418

Columns 12 through 22

-0.2130 -0.1878 -0.1655 -0.1453 -0.1269 -0.1100 -0.0941 -0.0792 -0.0650 -0.0514 -0.0382

Columns 23 through 25

-0.0253 -0.0126 0

>> idx=-N:N %harmonic indices

idx =

Columns 1 through 19

-25 -24 -23 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7

Columns 20 through 38

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12

Columns 39 through 51

13 14 15 16 17 18 19 20 21 22 23 24 25

>> Fna=5j./(idx*pi); %complex exponential terms

>> Fna(N+1)=5;

>> Bna=-2*imag(Fna(N+2:end)); %sine terms, actual ones (note that we didn’t calculate Ao because it’s neglibible (sawtooth is odd symmetric))

>> BnError=(Bn-Bna)./Bna %finding relative errors

BnError =

Columns 1 through 11

-0.0013 -0.0053 -0.0119 -0.0211 -0.0331 -0.0478 -0.0653 -0.0857 -0.1089 -0.1352 -0.1645

Columns 12 through 22

-0.1971 -0.2330 -0.2723 -0.3152 -0.3620 -0.4128 -0.4678 -0.5273 -0.5917 -0.6612 -0.7363

Columns 23 through 25

-0.8174 -0.9051 -1.0000

>> stem(idx,abs(Fn)) %plot line spectra of complex exponential Fourier series

>> axis tight

7. Control Flow

7.1 For loops

>> for n=1:10

x(n)=sin(n*pi/10);

end

>> x

x =

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000

>> for n=10:-1:1 %decrement

x(n)=sin(n*pi/10);

end

>> x

x =

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000

>> i=0;

>> for n=(1:10)’ %here right hand side of equal sign of n is a column, so the loop runs only once

i=i+1;

x(n)=sin(n*pi/10);

end

>> i

i =

>> x

x =

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000

>> array=randperm(10)

array =

6 3 7 8 5 1 2 4 9 10

>> for n=array %we can equate it to a variable array

x(n)=sin(n*pi/10);

end

>> x

x =

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000

>> for n=(1:10) %assigning n to 10 inside the loop won’t terminate it

x(n)=sin(n*pi/10);

n=10;

end

>> x

x =

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000

>> i=1

i =

>> r=rand(4,5)

r =

0.4173 0.4909 0.3692 0.2417 0.9421

0.0497 0.4893 0.1112 0.4039 0.9561

0.9027 0.3377 0.7803 0.0965 0.5752

0.9448 0.9001 0.3897 0.1320 0.0598

>> for x=r %here, y is an 1X5 array built as y(1) is sum of first column of x matrix (which is r), y(2) is sum of 2^nd column …

y(i)=sum(x);

i=i+1;

end

>> y

y =

2.3144 2.2179 1.6504 0.8740 2.5332

>> for n=1:5 %nested loops

for m=5:-1:1

A(n,m)=n^2+m^2;

end

disp(n)

end

>> A

A =

2 5 10 17 26

5 8 13 20 29

10 13 18 25 34

17 20 25 32 41

26 29 34 41 50

>> n=1:10;

>> x=sin(n*pi/10) %vectorized operation; much faster

x =

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000

>> n=1:5; m=1:5;

>> [nn,mm]=meshgrid(n,m);

>> A=nn.^2+mm.^2 %vectorized in 2D

A =

2 5 10 17 26

5 8 13 20 29

10 13 18 25 34

17 20 25 32 41

26 29 34 41 50

>> x=zeros(1,10); %preallocated memory for x; makes the loop much faster, otherwise it would increase size of x every iteration of the loop

>> for n=1:10

x(n)=sin(n*pi/10);

end

>> x

x =

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000

7.2 While loops

>> num=0; EPS=1;

>> while (EPS+1)>1

EPS=EPS/2;

num=num+1;

end

>> num

num =

>> EPS=EPS*2

EPS =

2.2204e-16

>> eps

ans =

2.2204e-16

7.3 If-Else-End Constructions

>> apples=10

apples =

>> cost=apples*25

cost =

250

>> if apples>5

cost=(1-20/100)*cost;

end

>> cost

cost =

200

>> EPS=1

EPS =

>> for num=1:100

EPS=EPS/2;

if(1+EPS)<=1

EPS=EPS*2

break

end

EPS =

2.2204e-16

>> num

num =

>> EPS=1

EPS =

>> for num=1:100

EPS=EPS/2;

if(1+EPS)>1

continue

end

EPS=EPS*2

break

end

EPS =

2.2204e-16

7.4 Switch-Case Constructions

>> x=2.7;

>> units=’m’;

>> switch units

case {‘inch’,’in’}

y=x*2.54;

case {‘feet’,’ft’}

y=x*2.54/12;

case {‘meter’,’m’}

y=x/100;

case {‘millimeter’,’mm’}

y=x*10;

case {‘centimeter’,’cm’}

y=x;

otherwise

disp([‘Unknown Units: ‘ units])

y=nan;

end

>> y

y =

0.0270

7.5 Try-Catch Blocks

>> x=ones(4,2)

x =

1 1

>> y=4*eye(2)

y =

4 0

0 4

>> try

z=x*y;

catch

z=nan;

disp(‘Error: X & Y are not conformable’)

end

>> z

z =

4 4

>> x=ones(4,2)

x =

1 1

>> y=4*eye(3)

y =

4 0 0

0 4 0

0 0 4

>> try

z=x*y;

catch

z=nan;

disp(‘Error: X & Y are not conformable’)

end

Error: X & Y are not conformable

>> lasterr

ans =

Error using *

Inner matrix dimensions must agree.

8. Functions

8.1 M-file Construction Rules

********************************M-FILE*************************

function d=mmempty(a,b)

%this is H1 line comes when lookfor command is used

%those lines are called when help command is executed

if isempty(a)

d=b;

else

d=a;

end

********************************M-FILE*************************

>> help mmempty

this is H1 line comes when lookfor command is used

those lines are called when help command is executed

>> a=[];

>> a

a =

[]

>> b=[1 4 3]

b =

1 4 3

>> mmempty(a,b)

ans =

1 4 3

>> a=[7 9 12]

a =

7 9 12

>> mmempty(a,b)

ans =

7 9 12

8.2 Input and Output Arguments

********************************M-FILE*************************

function y=mmdigit(x,n,b,t)

if nargin<2 %nargin is used to check number of inptut arguments; nargout for output arguments

error(‘Not enough input arguments’)

elseif nargin==2

b=10;

t=’round’;

elseif nargin==3

t=’round’;

end

n=round(abs(n(1)));

if isempty(b)

b=10;

else

b=round(abs(b(1)));

end

if isreal(x)

y=abs(x)+(x==0);

e=floor(log(y)./log(b)+1);

p=repmat(b,size(x)).^(n-e);

if strncmpi(t,’round’,1)

y=round(p.*x)./p;

elseif strncmpi(t,’fix’,2)

y=fix(p.*x)./p;

elseif strncmpi(t,’ceil’,1)

y=ceil(p.*x)./p;

elseif strncmpi(t,’floor’,2)

y=floor(p.*x)./p;

else

error(‘Unknown Rounding Requested’)

end

else

y=complex(mmdigit(real(x),n,b,t),mmdigit(imag(x),n,b,t))

end

********************************M-FILE*************************

function v=myfun(varargin)

%varargin function is used when we need unlimited number of input arguments; 1^st input goes to varargin{1} and second one goes to varagin{2} …

function v=myfun(x,y,z,varargin)

%here first input argument goes to x, second to y, third to z, fourth to varargin{1}, fifth to varargin{2} …

function v=myfun(varargout)

%varargout does the same for output arguments

********************************M-FILE*************************

8.3 Function Workspaces

********************************M-FILE*************************

function tic

global TICTOC %global makes it available for all function workspaces and MATLAB workspace

TICTOC=clock;

********************************M-FILE*************************

function t=toc

global TICTOC

if isempty(TICTOC)

error(‘Must call tic before toc’);

end

if nargout<1

elapsed_time=etime(clock,TICTOC)

else

t=etime(clock,TICTOC);

end

********************************M-FILE*************************

%persistent makes it available for functions workspace for all recursive calls for themselves

%evalin

%assignin

%inputname

8.4 Functions and the MATLAB search path

8.5 Creating your own toolbox

8.6 Command-Function Duality

8.7 Function evaluation using feval

>> fhan=@humps %using function handle to create fhan array as humps(x)

fhan =

@humps

>> fnc={@humps}

fnc =

@humps

>> fnc{2}=@cos

fnc =

@humps @cos

>> func2str(fhan) %returns function to a string

ans =

humps

>> feval(fhan,3) %evaluates the function at 3 ; note that feval works ONLY with M-file functions

ans =

-5.6383

>> feval(fnc{1},3)

ans =

-5.6383

>> feval(fnc{2},pi/2)

ans =

6.1232e-17

>> fncn={str2func(‘cos’)}

fncn =

@cos

>> fncn{2}=str2func(‘cos’)

fncn =

@cos @cos

>> fncn{1}=str2func(‘humps’)

fncn =

@humps @cos

>> fstr={‘humps’ ‘cos’}

fstr =

‘humps’ ‘cos’

>> funct=str2func(fstr)

Undefined function ‘str2func’ for input arguments of type ‘cell’.

>> funcc{1}=str2func(fstr{1})

funcc =

@humps

>> funcc{2}=str2func(fstr{2})

funcc =

@humps @cos

>> for(i=1:2)

fcn{i}=str2func(fstr{i})

end

fcn =

@humps @cos

>> functions(fhan) %gives information about fhan

ans =

function: ‘humps’

type: ‘simple’

file: ‘D:\MATLAB\toolbox\matlab\demos\humps.m’

>> myf=’100*(y-x^2)^2+(1-x)^2′ %create the function definition

myf =

100*(y-x^2)^2+(1-x)^2

>> x=3, y=5

x =

y =

>> a=eval(myf) %we can use eval to evaluate the function since it’s not an m-file function

a =

1604

>> b=feval(myf) %we can’t use feval here

Invalid function name ‘100*(y-x^2)^2+(1-x)^2’.

>> b=feval(myf,x,y) %even if you supply input arguments ; it’s not an m-file function

Invalid function name ‘100*(y-x^2)^2+(1-x)^2′.

>> myfi=inline(myf,’x’,’y’) %here’s one way to use feval to evalute it; inline function converts it to inline function

myfi =

Inline function:

myfi(x,y) = 100*(y-x^2)^2+(1-x)^2

>> b=feval(myfi,x,y) %and therefore be evaluated using feval

b =

1604

>> c=feval(myfi,-9,4) %even with different arguments

c =

593000

>> argnames(myfi) %check the arguments names

ans =

‘x’

‘y’

>> formula(myfi) %see the formula in myfi

ans =

100*(y-x^2)^2+(1-x)^2

>> fcnchk(fhan)

ans =

@humps

>> fun=’humps’

fun =

humps

>> fcnchk(fun)

ans =

@humps

>> fcnchk(myfi)

ans =

Inline function:

ans(x,y) = 100*(y-x^2)^2+(1-x)^2

>> fcnchk(myf)

ans =

Inline function:

ans(x,y) = 100*(y-x^2)^2+(1-x)^2

>> argnames(ans)

ans =

‘x’

‘y’

>> myfx=strrep(myf,’x’,’x(1)’)

myfx =

100*(y-x(1)^2)^2+(1-x(1))^2

>> myfy=strrep(myfx,’y’,’x(2)’)

myfy =

100*(x(2)-x(1)^2)^2+(1-x(1))^2

>> fcnchk(myfy)

ans =

Inline function:

ans(x) = 100*(x(2)-x(1)^2)^2+(1-x(1))^2

>> argnames(ans)

ans =

‘x’

9. Data Interpolation

9.1 1D Interpolation

Using the interp1 function (note that the independent variable MUST be monotonic; always decreasing or increasing)

>> x1=linspace(0,2*pi,60);

>> x2=linspace(0,2*pi,6);

>> plot(x1,sin(x1),x2,sin(x2),’–‘ %you can see the one with 60 points is more accurate than the one with 6 points

> Hz=[20:10:100 200:100:1000 1500 2000:1000:10000] %Frequencies in Hertz

Hz =

Columns 1 through 10

20 30 40 50 60 70 80 90 100 200

Columns 11 through 20

300 400 500 600 700 800 900 1000 1500 2000

Columns 21 through 28

3000 4000 5000 6000 7000 8000 9000 10000

>> spl=[76 66 59 54 49 46 43 40 38 22 14 9 6 3.5 2.5 1.4 0.7 0 -1 -3 -8 -7 -2 2 7 9 11 12] %sound pressure level in dB

spl =

Columns 1 through 11

76.0000 66.0000 59.0000 54.0000 49.0000 46.0000 43.0000 40.0000 38.0000 22.0000 14.0000

Columns 12 through 22

9.0000 6.0000 3.5000 2.5000 1.4000 0.7000 0 -1.0000 -3.0000 -8.0000 -7.0000

Columns 23 through 28

-2.0000 2.0000 7.0000 9.0000 11.0000 12.0000

>> plot(x1,sin(x1),x2,sin(x2),’–‘)

>> semilogx(Hz,spl,’-o’)

>> s=interp1(Hz,spl,2.5e3) %linear interpolation as default

s =

-5.5000

>> s=interp1(Hz,spl,2.5e3,’linear’) %linear interpolation; sufficient to most cases

s =

-5.5000

>> s=interp1(Hz,spl,2.5e3,’cubic’) %cubic interpolation

s =

-6.0488

>> s=interp1(Hz,spl,2.5e3,’spline’) %cubic spline; gives most accurate results but slowest method

s =

-5.8690

>> s=interp1(Hz,spl,2.5e3,’nearest’) %nearest neighbor method; fastest, but poor results

s =

-8

>> Hzi=linspace(2e3,5e3); %create data for x-axis near around the minimum value

>> spli=interp1(Hz,spl,Hzi,’spline’); %interpolate for the whole range of Hzi using cubic spline

>> i=find(Hz>=2e3 & Hz<=5e3) %find indices of original data points near minimum

i =

20 21 22 23

>> semilogx(Hz(i),spl(i),’–o’,Hzi,spli)

>> [spl_min,i]=min(spli) %find minimum spl value and its location

spl_min =

-8.4245

i =

>> Hz_min=Hzi(i) %find corresponding minimum Hz value

Hz_min =

3.3333e+03

>> x=linspace(0,2*pi,11)’

x =

0.6283

1.2566

1.8850

2.5133

3.1416

3.7699

4.3982

5.0265

5.6549

6.2832

>> y=[sin(x) cos(x) tan(x)]

y =

0 1.0000 0

0.5878 0.8090 0.7265

0.9511 0.3090 3.0777

0.9511 -0.3090 -3.0777

0.5878 -0.8090 -0.7265

0.0000 -1.0000 -0.0000

-0.5878 -0.8090 0.7265

-0.9511 -0.3090 3.0777

-0.9511 0.3090 -3.0777

-0.5878 0.8090 -0.7265

-0.0000 1.0000 -0.0000

>> xi=linspace(0,2*pi)

xi =

Columns 1 through 11

0 0.0635 0.1269 0.1904 0.2539 0.3173 0.3808 0.4443 0.5077 0.5712 0.6347

Columns 12 through 22

0.6981 0.7616 0.8251 0.8885 0.9520 1.0155 1.0789 1.1424 1.2059 1.2693 1.3328

Columns 23 through 33

1.3963 1.4597 1.5232 1.5867 1.6501 1.7136 1.7771 1.8405 1.9040 1.9675 2.0309

Columns 34 through 44

2.0944 2.1579 2.2213 2.2848 2.3483 2.4117 2.4752 2.5387 2.6021 2.6656 2.7291

Columns 45 through 55

2.7925 2.8560 2.9195 2.9829 3.0464 3.1099 3.1733 3.2368 3.3003 3.3637 3.4272

Columns 56 through 66

3.4907 3.5541 3.6176 3.6811 3.7445 3.8080 3.8715 3.9349 3.9984 4.0619 4.1253

Columns 67 through 77

4.1888 4.2523 4.3157 4.3792 4.4427 4.5061 4.5696 4.6331 4.6965 4.7600 4.8235

Columns 78 through 88

4.8869 4.9504 5.0139 5.0773 5.1408 5.2043 5.2677 5.3312 5.3947 5.4581 5.5216

Columns 89 through 99

5.5851 5.6485 5.7120 5.7755 5.8389 5.9024 5.9659 6.0293 6.0928 6.1563 6.2197

Column 100

6.2832

>> yi=interp1(x,y,xi,’cubic’) %interpolate on a finer data interval

yi =

0 1.0000 0

0.0698 0.9941 0.0106

0.1377 0.9841 0.0408

0.2035 0.9703 0.0887

0.2671 0.9533 0.1520

0.3283 0.9334 0.2286

0.3869 0.9111 0.3165

0.4427 0.8868 0.4135

0.4957 0.8610 0.5174

0.5456 0.8340 0.6262

0.5923 0.8062 0.7382

0.6397 0.7737 0.9046

0.6890 0.7343 1.1452

0.7387 0.6890 1.4378

0.7869 0.6390 1.7601

0.8321 0.5853 2.0900

0.8723 0.5289 2.4052

0.9061 0.4709 2.6835

0.9315 0.4124 2.9028

0.9470 0.3544 3.0408

0.9511 0.2978 3.0702

0.9511 0.2394 2.8283

0.9511 0.1779 2.3009

0.9511 0.1143 1.5641

0.9511 0.0492 0.6941

0.9511 -0.0164 -0.2330

0.9511 -0.0819 -1.1410

0.9511 -0.1463 -1.9539

0.9511 -0.2090 -2.5955

0.9511 -0.2690 -2.9897

0.9505 -0.3259 -3.0723

0.9406 -0.3833 -2.9833

0.9200 -0.4417 -2.8019

0.8901 -0.5001 -2.5504

0.8529 -0.5574 -2.2508

0.8100 -0.6126 -1.9255

0.7631 -0.6646 -1.5966

0.7139 -0.7123 -1.2864

0.6642 -0.7548 -1.0170

0.6156 -0.7909 -0.8108

0.5692 -0.8202 -0.6829

0.5191 -0.8482 -0.5838

0.4646 -0.8757 -0.4971

0.4065 -0.9020 -0.4202

0.3456 -0.9266 -0.3510

0.2828 -0.9486 -0.2870

0.2190 -0.9676 -0.2259

0.1550 -0.9828 -0.1652

0.0918 -0.9936 -0.1026

0.0300 -0.9993 -0.0358

-0.0300 -0.9993 0.0358

-0.0918 -0.9936 0.1026

-0.1550 -0.9828 0.1652

-0.2190 -0.9676 0.2259

-0.2828 -0.9486 0.2870

-0.3456 -0.9266 0.3510

-0.4065 -0.9020 0.4202

-0.4646 -0.8757 0.4971

-0.5191 -0.8482 0.5838

-0.5692 -0.8202 0.6829

-0.6156 -0.7909 0.8108

-0.6642 -0.7548 1.0170

-0.7139 -0.7123 1.2864

-0.7631 -0.6646 1.5966

-0.8100 -0.6126 1.9255

-0.8529 -0.5574 2.2508

-0.8901 -0.5001 2.5504

-0.9200 -0.4417 2.8019

-0.9406 -0.3833 2.9833

-0.9505 -0.3259 3.0723

-0.9511 -0.2690 2.9897

-0.9511 -0.2090 2.5955

-0.9511 -0.1463 1.9539

-0.9511 -0.0819 1.1410

-0.9511 -0.0164 0.2330

-0.9511 0.0492 -0.6941

-0.9511 0.1143 -1.5641

-0.9511 0.1779 -2.3009

-0.9511 0.2394 -2.8283

-0.9511 0.2978 -3.0702

-0.9470 0.3544 -3.0408

-0.9315 0.4124 -2.9028

-0.9061 0.4709 -2.6835

-0.8723 0.5289 -2.4052

-0.8321 0.5853 -2.0900

-0.7869 0.6390 -1.7601

-0.7387 0.6890 -1.4378

-0.6890 0.7343 -1.1452

-0.6397 0.7737 -0.9046

-0.5923 0.8062 -0.7382

-0.5456 0.8340 -0.6262

-0.4957 0.8610 -0.5174

-0.4427 0.8868 -0.4135

-0.3869 0.9111 -0.3165

-0.3283 0.9334 -0.2286

-0.2671 0.9533 -0.1520

-0.2035 0.9703 -0.0887

-0.1377 0.9841 -0.0408

-0.0698 0.9941 -0.0106

-0.0000 1.0000 -0.0000

9.2 2D Interpolation

Using the interp2 function

>> x=0:0.5:4

x =

0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000

>> y=0:0.5:6

y =

Columns 1 through 12

0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 4.5000 5.0000 5.5000

Column 13

6.0000

>> z=[100 99 100 99 100 99 99 99 100

100 99 99 99 100 99 100 99 99

99 99 98 98 100 99 100 100 100

100 98 97 97 99 100 100 100 99

101 100 98 98 100 102 103 100 100

102 103 101 100 102 106 104 101 100

99 102 100 100 103 108 106 101 99

97 99 100 100 102 105 103 101 100

100 102 103 101 102 103 102 100 99

100 102 103 102 101 101 100 99 99

100 100 101 101 100 100 100 99 99

100 100 100 100 100 99 99 99 99

100 100 100 99 99 100 99 100 99]

z =

100 99 100 99 100 99 99 99 100

100 99 99 99 100 99 100 99 99

99 99 98 98 100 99 100 100 100

100 98 97 97 99 100 100 100 99

101 100 98 98 100 102 103 100 100

102 103 101 100 102 106 104 101 100

99 102 100 100 103 108 106 101 99

97 99 100 100 102 105 103 101 100

100 102 103 101 102 103 102 100 99

100 102 103 102 101 101 100 99 99

100 100 101 101 100 100 100 99 99

100 100 100 100 100 99 99 99 99

100 100 100 99 99 100 99 100 99

>> mesh(x,y,z)

>> zi=interp2(x,y,z,2.2,3.3)

zi =

103.9200

>> zi=interp2(x,y,z,2.2,3.3,’linear’)

zi =

103.9200

>> zi=interp2(x,y,z,2.2,3.3,’cubic’)

zi =

104.1861

>> zi=interp2(x,y,z,2.2,3.3,’nearest’)

zi =

102

>> xi=linspace(0,4,30);

>> yi=linspace(0,6,40);

>> [xxi,yyi]=meshgrid(xi,yi); %to have all combinations of xi and yi, see below more illustration

>> zzi=interp2(x,y,z,xxi,yyi,’cubic’); %interpolate over an interval

>> mesh(xxi,yyi,zzi)

>> mesh(xxi,yyi,zzi), hold on, [xx,yy]=meshgrid(x,y), plot3(xx,yy,z+0.1,’ok’) %or can do this to show data nodes

>> zmax=max(max(zzi)) %fniding z maximum

zmax =

108.0520

>> [i,j]=find(zmax==zzi) %finding location of z maximum

i =

j =

>> xmax=xi(i)

xmax =

2.6207

>> ymax=yi(i)

ymax =

2.9231

>> xtest=1:5, ytest= 6:9, [xx,yy]=meshgrid(xtest,ytest) %illustrating meshgrid: it creates xx by creating number of rows = length(ytest) and yy by transposing y and creating number of columns = length(xtest)

xtest =

1 2 3 4 5

ytest =

6 7 8 9

xx =

1 2 3 4 5

yy =

6 6 6 6 6

7 7 7 7 7

8 8 8 8 8

9 9 9 9 9

For higher order matrices, use the interp3 and interpn functions and ndgrid instead of meshgrid

9.3 Triangulation and Scattered Data

Using the delaunay, trimesh, delaunayTriangulation, pointLocation, dsearchn, convhull, voronoi, griddata functions

>> x=randn(1,12)

x =

Columns 1 through 11

0.5377 1.8339 -2.2588 0.8622 0.3188 -1.3077 -0.4336 0.3426 3.5784 2.7694 -1.3499

Column 12

3.0349

>> y=randn(1,12)

y =

Columns 1 through 11

0.7254 -0.0631 0.7147 -0.2050 -0.1241 1.4897 1.4090 1.4172 0.6715 -1.2075 0.7172

Column 12

1.6302

>> plot(x,y,’o’)

>> tri=delaunay(x,y)

tri =

7 6 11

11 3 5

3 11 6

7 11 5

1 7 5

8 6 7

1 8 7

12 6 8

4 2 1

5 4 1

1 12 8

1 2 12

3 10 5

2 9 12

10 2 4

10 9 2

5 10 4

>> hold on, trimesh(tri,x,y)

>> DT = delaunayTriangulation(x’,y’)

DT =

delaunayTriangulation with properties:

Points: [12×2 double]

ConnectivityList: [17×3 double]

Constraints: []

>> pointLocation(DT,0,0)

ans =

>> DT(ans,:)

ans =

11 3 5

>> pointLocation(DT,[0,0;0,-0.5])

ans =

NaN

>> dsearchn([x’ y’],tri,[0 0.5; 1 2])

ans =

>> k=convhull(x,y)’

k =

3 10 9 12 6 3

>> hold on, plot(x,y,’o’,x(k),y(k))

>> voronoi(x,y,tri)

>> z=rand(1,12)

z =

Columns 1 through 11

0.6787 0.7577 0.7431 0.3922 0.6555 0.1712 0.7060 0.0318 0.2769 0.0462 0.0971

Column 12

0.8235

>> xi=linspace(min(x),max(x),30)

xi =

Columns 1 through 11

-2.2588 -2.0576 -1.8563 -1.6550 -1.4537 -1.2524 -1.0511 -0.8499 -0.6486 -0.4473 -0.2460

Columns 12 through 22

-0.0447 0.1566 0.3578 0.5591 0.7604 0.9617 1.1630 1.3643 1.5656 1.7668 1.9681

Columns 23 through 30

2.1694 2.3707 2.5720 2.7733 2.9745 3.1758 3.3771 3.5784

>> yi=linspace(min(y),max(y),30)

yi =

Columns 1 through 11

-1.2075 -1.1096 -1.0118 -0.9139 -0.8161 -0.7182 -0.6204 -0.5225 -0.4247 -0.3268 -0.2290

Columns 12 through 22

-0.1311 -0.0333 0.0646 0.1624 0.2603 0.3582 0.4560 0.5539 0.6517 0.7496 0.8474

Columns 23 through 30

0.9453 1.0431 1.1410 1.2388 1.3367 1.4345 1.5324 1.6302

>> [Xi,Yi]=meshgrid(xi,yi);

>> Zi=griddata(x,y,z,Xi,Yi);

>> mesh(Xi,Yi,Zi)

>> hold on, plot3(x,y,z,’ko’)

10. Optimization

10.1 Zero Finding

>> x=linspace(-.5,1.5);

>> y=humps(x);

>> plot(x,y), grid on

>> format long

>> x=fzero(‘humps’,1.3) %fzero looks in the function for x where f=0 given the function name and initial guess

x =

1.299549682584822

>> [x,val]=fzero(‘humps’,-0.2) %returns x where f=0 and the value at f to check acuuracy

x =

-0.131618018099606

val =

8.881784197001252e-16

>> [x,val]=fzero(‘humps’,[-2 0]) %we can provide fzero function with an interval, not just a single number

x =

-0.131618018099606

val =

>> [x,val]=fzero(‘humps’,[0 1.2]) %here we got an error because there’s no zero in this interval

Error using fzero (line 274)

The function values at the interval endpoints must differ in sign.

>> [x,val]=fzero(@humps,[-2 0]) %since humps is an m-file, we can use function handle approach

x =

-0.131618018099606

val =

>> humpstr=’1./((x-.3)^2+.01)+1./((x-.9)^2+.04)-6′

humpstr =

1./((x-.3)^2+.01)+1./((x-.9)^2+.04)-6

>> [x,val]=fzero(humpstr,[-2 0]) %string exprssion works also

x =

-0.131618018099606

val =

>> hinline=inline(humpstr)

hinline =

Inline function:

hinline(x) = 1./((x-.3)^2+.01)+1./((x-.9)^2+.04)-6

>> [x,val]=fzero(hinline,[-2 0]) %inline function object also works; exactly same as above method; above method: fzero makes the conversion to in-line function object, contrary here, we do it ourselves. NOTE: all approaches give same results, but function handle is fastest, then using function m-file directly, then either inline object or string expression

x =

-0.131618018099606

val =

>> options=optimset(‘Display’,’iter’) %show iteration history

options =

Display: ‘iter’

MaxFunEvals: []

MaxIter: []

TolFun: []

TolX: []

FunValCheck: []

OutputFcn: []

PlotFcns: []

ActiveConstrTol: []

Algorithm: []

AlwaysHonorConstraints: []

BranchStrategy: []

DerivativeCheck: []

Diagnostics: []

DiffMaxChange: []

DiffMinChange: []

FinDiffRelStep: []

FinDiffType: []

GoalsExactAchieve: []

GradConstr: []

GradObj: []

HessFcn: []

Hessian: []

HessMult: []

HessPattern: []

HessUpdate: []

InitialHessType: []

InitialHessMatrix: []

InitBarrierParam: []

InitTrustRegionRadius: []

Jacobian: []

JacobMult: []

JacobPattern: []

LargeScale: []

MaxNodes: []

MaxPCGIter: []

MaxProjCGIter: []

MaxRLPIter: []

MaxSQPIter: []

MaxTime: []

MeritFunction: []

MinAbsMax: []

NodeDisplayInterval: []

NodeSearchStrategy: []

NoStopIfFlatInfeas: []

ObjectiveLimit: []

PhaseOneTotalScaling: []

Preconditioner: []

PrecondBandWidth: []

RelLineSrchBnd: []

RelLineSrchBndDuration: []

ScaleProblem: []

Simplex: []

SubproblemAlgorithm: []

TolCon: []

TolConSQP: []

TolGradCon: []

TolPCG: []

TolProjCG: []

TolProjCGAbs: []

TolRLPFun: []

TolXInteger: []

TypicalX: []

UseParallel: []

>> [x,val]=fzero(@humps,[-2 0],options) %let’s see interpolation process!

Func-count x f(x) Procedure

2 0 5.17647 initial

3 -0.952481 -5.07853 interpolation

4 -0.480789 -3.87242 interpolation

5 -0.240394 -1.94304 bisection

6 -0.120197 0.28528 bisection

7 -0.135585 -0.0944316 interpolation

8 -0.131759 -0.00338409 interpolation

9 -0.131618 1.63632e-06 interpolation

10 -0.131618 -7.14819e-10 interpolation

11 -0.131618 0 interpolation

Zero found in the interval [-2, 0]

x =

-0.131618018099606

val =

>> options=optimset(‘Display’,’none’); %show nothing

>> [x,val]=fzero(@humps,[-2 0],options)

x =

-0.131618018099606

val =

>> options=optimset(‘Display’,’iter’,’TolX’,0.1); %show iteration history and set tolerance to 0.1

>> [x,val]=fzero(@humps,[-2 0],options)

Func-count x f(x) Procedure

2 0 5.17647 initial

3 -0.952481 -5.07853 interpolation

4 -0.480789 -3.87242 interpolation

5 -0.240394 -1.94304 bisection

Zero found in the interval [-2, 0]

x =

-0.240394472507622

val =

-1.943038259725652

>> help optimset

optimset Create/alter optimization OPTIONS structure.

OPTIONS = optimset(‘PARAM1′,VALUE1,’PARAM2′,VALUE2,…) creates an

optimization options structure OPTIONS in which the named parameters have

the specified values. Any unspecified parameters are set to [] (parameters

with value [] indicate to use the default value for that parameter when

OPTIONS is passed to the optimization function). It is sufficient to type

only the leading characters that uniquely identify the parameter. Case is

ignored for parameter names.

NOTE: For values that are strings, the complete string is required.

OPTIONS = optimset(OLDOPTS,’PARAM1’,VALUE1,…) creates a copy of OLDOPTS

with the named parameters altered with the specified values.

OPTIONS = optimset(OLDOPTS,NEWOPTS) combines an existing options structure

OLDOPTS with a new options structure NEWOPTS. Any parameters in NEWOPTS

with non-empty values overwrite the corresponding old parameters in

OLDOPTS.

optimset with no input arguments and no output arguments displays all

parameter names and their possible values, with defaults shown in {}

when the default is the same for all functions that use that parameter.

Use optimset(OPTIMFUNCTION) to see parameters for a specific function.

OPTIONS = optimset (with no input arguments) creates an options structure

OPTIONS where all the fields are set to [].

OPTIONS = optimset(OPTIMFUNCTION) creates an options structure with all

the parameter names and default values relevant to the optimization

function named in OPTIMFUNCTION. For example,

optimset(‘fminbnd’)

optimset(@fminbnd)

returns an options structure containing all the parameter names and

default values relevant to the function ‘fminbnd’.

optimset PARAMETERS for MATLAB

Display – Level of display [ off | iter | notify | final ]

MaxFunEvals – Maximum number of function evaluations allowed

[ positive integer ]

MaxIter – Maximum number of iterations allowed [ positive scalar ]

TolFun – Termination tolerance on the function value [ positive scalar ]

TolX – Termination tolerance on X [ positive scalar ]

FunValCheck – Check for invalid values, such as NaN or complex, from

user-supplied functions [ {off} | on ]

OutputFcn – Name(s) of output function [ {[]} | function ]

All output functions are called by the solver after each

iteration.

PlotFcns – Name(s) of plot function [ {[]} | function ]

Function(s) used to plot various quantities in every iteration

Note to Optimization Toolbox users:

To see the parameters for a specific function, check the documentation page

for that function. For instance, enter

doc fmincon

to open the reference page for fmincon.

You can also see the options in the Optimization Tool. Enter

optimtool

Examples:

To create an options structure with the default parameters for FZERO

options = optimset(‘fzero’);

To create an options structure with TolFun equal to 1e-3

options = optimset(‘TolFun’,1e-3);

To change the Display value of options to ‘iter’

options = optimset(options,’Display’,’iter’);

See also optimget, fzero, fminbnd, fminsearch, lsqnonneg.

Overloaded methods:

cgoptimstore/optimset

ResponseOptimizer.optimset

ParameterEstimator.optimset

10.2 Minimization in one dimension

>> [xmin,value]=fminbnd(‘humps’,0.5,0.8) %fminbnd finds the minimum value in an interval (here between 0.5 and 0.8)

xmin =

0.637008211963619

value =

11.252754125877694

>> options=optimset(‘Display’,’iter’)

options =

Display: ‘iter’

MaxFunEvals: []

MaxIter: []

TolFun: []

TolX: []

FunValCheck: []

OutputFcn: []

PlotFcns: []

ActiveConstrTol: []

Algorithm: []

AlwaysHonorConstraints: []

BranchStrategy: []

DerivativeCheck: []

Diagnostics: []

DiffMaxChange: []

DiffMinChange: []

FinDiffRelStep: []

FinDiffType: []

GoalsExactAchieve: []

GradConstr: []

GradObj: []

HessFcn: []

Hessian: []

HessMult: []

HessPattern: []

HessUpdate: []

InitialHessType: []

InitialHessMatrix: []

InitBarrierParam: []

InitTrustRegionRadius: []

Jacobian: []

JacobMult: []

JacobPattern: []

LargeScale: []

MaxNodes: []

MaxPCGIter: []

MaxProjCGIter: []

MaxRLPIter: []

MaxSQPIter: []

MaxTime: []

MeritFunction: []

MinAbsMax: []

NodeDisplayInterval: []

NodeSearchStrategy: []

NoStopIfFlatInfeas: []

ObjectiveLimit: []

PhaseOneTotalScaling: []

Preconditioner: []

PrecondBandWidth: []

RelLineSrchBnd: []

RelLineSrchBndDuration: []

ScaleProblem: []

Simplex: []

SubproblemAlgorithm: []

TolCon: []

TolConSQP: []

TolGradCon: []

TolPCG: []

TolProjCG: []

TolProjCGAbs: []

TolRLPFun: []

TolXInteger: []

TypicalX: []

UseParallel: []

>> [xmin,value]=fminbnd(‘humps’,0.5,0.8,options) %let’s see the iterations

Func-count x f(x) Procedure

1 0.61459 11.4103 initial

2 0.68541 11.9288 golden

3 0.57082 12.7389 golden

4 0.638866 11.2538 parabolic

5 0.637626 11.2529 parabolic

6 0.637046 11.2528 parabolic

7 0.637008 11.2528 parabolic

8 0.636975 11.2528 parabolic

Optimization terminated:

the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04

xmin =

0.637008211963619

value =

11.252754125877694

>> [xmax,value]=fminbnd(‘-humps(x)’,0.2,0.4,options) %let’s find the maximum using the string expression approach, or alternatively we could negate the expression in the humps m-file, note that the real maximum value is the negative of the value given here

Func-count x f(x) Procedure

1 0.276393 -91.053 initial

2 0.323607 -91.4079 golden

3 0.352786 -75.1541 golden

4 0.300509 -96.5012 parabolic

5 0.300397 -96.5014 parabolic

6 0.300364 -96.5014 parabolic

7 0.300331 -96.5014 parabolic

Optimization terminated:

the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04

xmax =

0.300364137900245

value =

-96.501407243870503

10.3 Minimization in higher dimensions

>> x=[-1.5:0.125:1.5];

>> y=[-.6:0.125:2.8];

>> [X,Y]=meshgrid(x,y);

>> Z=100.*(Y-X.*X).^2+(1-X).^2;

>> mesh(X,Y,Z)

>> x=[-1.5:0.125:1.5];

>> y=[-.6:0.125:2.8];

>> [X,Y]=meshgrid(x,y);

>> Z=100.*(Y-X.*X).^2+(1-X).^2;

>> mesh(X,Y,Z)

********************************M-FILE*************************

function f=banana(x)

f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;

********************************M-FILE*************************

>> [xmin,value,flag,output]=fminsearch(‘banana’,[-1.9,2])

xmin =

1.000016668894802 1.000034473862770

value =

4.068551535063419e-10

flag =

output =

iterations: 114

funcCount: 210

algorithm: ‘Nelder-Mead simplex direct search’

message: [1×194 char]

>> options=optimset(‘Display’,’none’,’TolFun’,1e-8,’TolX’,1e-8)

options =

Display: ‘none’

MaxFunEvals: []

MaxIter: []

TolFun: 1.000000000000000e-08

TolX: 1.000000000000000e-08

FunValCheck: []

OutputFcn: []

PlotFcns: []

ActiveConstrTol: []

Algorithm: []

AlwaysHonorConstraints: []

BranchStrategy: []

DerivativeCheck: []

Diagnostics: []

DiffMaxChange: []

DiffMinChange: []

FinDiffRelStep: []

FinDiffType: []

GoalsExactAchieve: []

GradConstr: []

GradObj: []

HessFcn: []

Hessian: []

HessMult: []

HessPattern: []

HessUpdate: []

InitialHessType: []

InitialHessMatrix: []

InitBarrierParam: []

InitTrustRegionRadius: []

Jacobian: []

JacobMult: []

JacobPattern: []

LargeScale: []

MaxNodes: []

MaxPCGIter: []

MaxProjCGIter: []

MaxRLPIter: []

MaxSQPIter: []

MaxTime: []

MeritFunction: []

MinAbsMax: []

NodeDisplayInterval: []

NodeSearchStrategy: []

NoStopIfFlatInfeas: []

ObjectiveLimit: []

PhaseOneTotalScaling: []

Preconditioner: []

PrecondBandWidth: []

RelLineSrchBnd: []

RelLineSrchBndDuration: []

ScaleProblem: []

Simplex: []

SubproblemAlgorithm: []

TolCon: []

TolConSQP: []

TolGradCon: []

TolPCG: []

TolProjCG: []

TolProjCGAbs: []

TolRLPFun: []

TolXInteger: []

TypicalX: []

UseParallel: []

>> [xmin,value,flag,output]=fminsearch(‘banana’,[-1.9,2],options)

xmin =

1.000000001260770 1.000000002307897

value =

6.153858843361103e-18

flag =

output =

iterations: 144

funcCount: 266

algorithm: ‘Nelder-Mead simplex direct search’

message: [1×194 char]

10.4 Practical Issues

11. 2D Graphics

11.1 The plot function

>> x=linspace(0,2*pi,30),y=sin(x), z=cos(x)

x =

Columns 1 through 12

0 0.2167 0.4333 0.6500 0.8666 1.0833 1.3000 1.5166 1.7333 1.9500 2.1666 2.3833

Columns 13 through 24

2.5999 2.8166 3.0333 3.2499 3.4666 3.6832 3.8999 4.1166 4.3332 4.5499 4.7666 4.9832

Columns 25 through 30

5.1999 5.4165 5.6332 5.8499 6.0665 6.2832

y =

Columns 1 through 12

0 0.2150 0.4199 0.6052 0.7622 0.8835 0.9635 0.9985 0.9868 0.9290 0.8277 0.6877

Columns 13 through 24

0.5156 0.3193 0.1081 -0.1081 -0.3193 -0.5156 -0.6877 -0.8277 -0.9290 -0.9868 -0.9985 -0.9635

Columns 25 through 30

-0.8835 -0.7622 -0.6052 -0.4199 -0.2150 -0.0000

z =

Columns 1 through 12

1.0000 0.9766 0.9076 0.7961 0.6474 0.4684 0.2675 0.0541 -0.1618 -0.3701 -0.5612 -0.7260

Columns 13 through 24

-0.8569 -0.9477 -0.9941 -0.9941 -0.9477 -0.8569 -0.7260 -0.5612 -0.3701 -0.1618 0.0541 0.2675

Columns 25 through 30

0.4684 0.6474 0.7961 0.9076 0.9766 1.0000

>> plot(x,y,x,z), title(‘Sine and Cosine’) %plotting and customizing the title

>> W=[y; z]

W =

Columns 1 through 11

0 0.2150 0.4199 0.6052 0.7622 0.8835 0.9635 0.9985 0.9868 0.9290 0.8277

1.0000 0.9766 0.9076 0.7961 0.6474 0.4684 0.2675 0.0541 -0.1618 -0.3701 -0.5612

Columns 12 through 22

0.6877 0.5156 0.3193 0.1081 -0.1081 -0.3193 -0.5156 -0.6877 -0.8277 -0.9290 -0.9868

-0.7260 -0.8569 -0.9477 -0.9941 -0.9941 -0.9477 -0.8569 -0.7260 -0.5612 -0.3701 -0.1618

Columns 23 through 30

-0.9985 -0.9635 -0.8835 -0.7622 -0.6052 -0.4199 -0.2150 -0.0000

0.0541 0.2675 0.4684 0.6474 0.7961 0.9076 0.9766 1.0000

>> plot(x,W) %same as before

>> plot(W,x), title(‘Changing arguments order’)

>> plot(x) %plot can take 1 argument; it plots index of x versus value

>> d=randn(5,1)+i*randn(5,1) %building a complex numbers vector

d =

-1.1471 + 0.3252i

-1.0689 – 0.7549i

-0.8095 + 1.3703i

-2.9443 – 1.7115i

1.4384 – 0.1022i

>> plot(d) %here MATLAB plots real(d) vs imag(d)

11.2 Linestyles, Markers, Colors

>> plot(x,y,’b:p’,x,z,’c-‘) %change color, linestyle, marker : color: b;blue g;green r:red c:cyan m:magenta y:yellow k:black w:white linestyle: . o x + * s:square d:diamond v:triangle up ^triangle down <triangle right >triangle left p pentagram hhexagram markers: – : -. —

11.3 Plot Grids, Axes Box, Labels

>> plot(x,y,x,z), box off, xlabel(‘X’),ylabel(‘Y and Z’), title(‘Sine & Cosine’) %turning box off, adding labels to axes

>> plot(x,y,x,z), grid on, text(2.5,0.7,’sin(x)’) %adding grid lines, adding text at (x,y) of (2.5,7) as “sin(X)”

11.4 Customizing Plot Axes

>> plot(x,y), axis([0 2*pi -1.5 2]) %control the limits of axes, if you wish to set only x or y we can use xlim or ylim ([p1 p2]) note that there are many axis options and settings: V=axis return axis limits, axis auto: makes it auto, axis manual: freeze axis scaling (if hold on is used, next plots forced to use such scale), axis tight: set axis limits to range of plotted data, axis fill: set axis limits and aspect ratio so that axis fills allotted space(only if plotboxaspectraio or dataaspectratiomode is on), axis ij: axes in matrix mode; X increase from left to right and y from top to bottom, axis xy axes in Cartesian mode; X increase from left to right and y from bottom to top, axis equal, axis image, axis square, axis normal, axis vis3d: freeze aspect ratio to enable 3D object rotation, axis off, axis on

11.5 Multiple Plots

>> plot(x,y), hold on

>> ishold

ans =

>> plot(x,z,’m’)

>> hold off

>> ishold

ans =

>> title (‘playing with hold’)

11.6 Multiple Figures

>> plot(x,y), figure, plot(x,z) %plots (x,z) in a new figure

>> close(2) %close figure(2)

>> close all %close all

>> plot(x,y), figure, plot(x,z)

>> clf //clear current figure

>> clf reset //clear current figure and reset all properties(like hold) to their default values

11.7 Subplots

>> x=linspace(0,2*pi,30), y=sin(x), z=cos(x), a=2*sin(x).*cos(x), b=sin(x)./(cos(x)+eps)

x =

Columns 1 through 12

0 0.2167 0.4333 0.6500 0.8666 1.0833 1.3000 1.5166 1.7333 1.9500 2.1666 2.3833

Columns 13 through 24

2.5999 2.8166 3.0333 3.2499 3.4666 3.6832 3.8999 4.1166 4.3332 4.5499 4.7666 4.9832

Columns 25 through 30

5.1999 5.4165 5.6332 5.8499 6.0665 6.2832

y =

Columns 1 through 12

0 0.2150 0.4199 0.6052 0.7622 0.8835 0.9635 0.9985 0.9868 0.9290 0.8277 0.6877

Columns 13 through 24

0.5156 0.3193 0.1081 -0.1081 -0.3193 -0.5156 -0.6877 -0.8277 -0.9290 -0.9868 -0.9985 -0.9635

Columns 25 through 30

-0.8835 -0.7622 -0.6052 -0.4199 -0.2150 -0.0000

z =

Columns 1 through 12

1.0000 0.9766 0.9076 0.7961 0.6474 0.4684 0.2675 0.0541 -0.1618 -0.3701 -0.5612 -0.7260

Columns 13 through 24

-0.8569 -0.9477 -0.9941 -0.9941 -0.9477 -0.8569 -0.7260 -0.5612 -0.3701 -0.1618 0.0541 0.2675

Columns 25 through 30

0.4684 0.6474 0.7961 0.9076 0.9766 1.0000

a =

Columns 1 through 12

0 0.4199 0.7622 0.9635 0.9868 0.8277 0.5156 0.1081 -0.3193 -0.6877 -0.9290 -0.9985

Columns 13 through 24

-0.8835 -0.6052 -0.2150 0.2150 0.6052 0.8835 0.9985 0.9290 0.6877 0.3193 -0.1081 -0.5156

Columns 25 through 30

-0.8277 -0.9868 -0.9635 -0.7622 -0.4199 -0.0000

b =

Columns 1 through 12

0 0.2201 0.4626 0.7602 1.1773 1.8862 3.6017 18.4439 -6.0997 -2.5098 -1.4749 -0.9473

Columns 13 through 24

-0.6017 -0.3369 -0.1088 0.1088 0.3369 0.6017 0.9473 1.4749 2.5098 6.0997 -18.4439 -3.6017

Columns 25 through 30

-1.8862 -1.1773 -0.7602 -0.4626 -0.2201 -0.0000

>> subplot(2,2,1), plot(x,y), axis([0 2*pi -1 1]) %subplot tells MATLAB to divide the figure into 2 X 2 matrix, to activate a specific part, you type 1 to select the part on the left on the first row; MATLAB orders them by left to right in top row, then next rown and so on, note that all graphics options affect only current subplot; to return back to 1 plot in the figure, type subplot(1,1,1)

>> subplot(2,2,2), plot(x,z), axis([0 2*pi -1 1])

>> subplot(2,2,3), plot(x,a), axis([0 2*pi -1 1])

>> subplot(2,2,4), plot(x,b), axis([0 2*pi -20 20])

11.8 Interactive Plotting Tools

>> plot(x,y,x,z), legend(‘sin(x)’,’cos(x)’) %adding a legend

>> plot(x,y,x,z), legend(‘sin(x)’,’cos(x)’), zoom(3) %causes the plot to be zoomed in 3 times

>> plot(x,y,x,z), legend(‘sin(x)’,’cos(x)’), [l,m]=ginput(5) %takes location(x,y) of 5 points clicked on plot

l =

1.3952

-1.1532

1.1532

2.7661

m =

0.1842

0.9152

-0.1140

0.3830

11.9 Screen Updates

>> plot(x,y,x,z), legend(‘sin(x)’,’cos(x)’), grid %writing legend() and grid in same line forces MATLAB to update the screenl if written on separate lines, it would not update until next plot command; 5 events can update a figure: enter is pressed, function that stops execution (pause, keyboard, input, …), getframe, drawnow (only way to force MATLAB to update screen), resize figure window

11.10 Specialized 2D Plots

>> c=-pi:pi/5:pi

c =

-3.1416 -2.5133 -1.8850 -1.2566 -0.6283 0 0.6283 1.2566 1.8850 2.5133 3.1416

>> area([sin(c);cos(c)]) %area() plots a column vector and fills the area between the curve and X-axis, if it’s a matrix, like the situation here, it plots each column individually and again fills area between each curve and X-axis (note that it adds values of previous column when plotting next column; for instance here, first curve is only 2 points sin(-pi) and cos(-pi) which is (0,1); it plots 0 at x = 1 and 1 at x = 2 and second curve is at sin(-2.5133) and cos(-2.5133) which is (-0.5877,-0.809) adding to previous point (0,-1) yields (-0.58,-1.80) which is overlaid by next curves, third curve happens at third column sin(-1.885) and cos(-1.885) which is (-.95104,-.30905) adding to previous columns (0-.58-.95,-1-.80-0.30) which is (-1.53,-2.1) which has an overlaid part and you can see the blue part; fourth happens at (0-0.5877-0.95104-0.95104,-1-0.80903-0.30905+0.30905) or (-2.48,-1.81) which is overlaid; fifth happens at (-2.48-0.5877,-1.81+0.81) or (-3.06,-1) so the curve is -3.06 at x = 0 and -1 at x = 2 that’s the green area; sixth curve is at (-3.06+0,-1+1) or (-3.06,0) which is the yellow one; seventh curve is at (-3.06+0.5877,0+0.809) or (-2.47,0.81) which is the orange curve; eighth curve is at (-2.47+0.95,0.81+.309) or (-1.52,1.11) which is the orange one

>> t=(1:2:15)’*pi/8,x=sin(t),y=cos(t)

t =

0.3927

1.1781

1.9635

2.7489

3.5343

4.3197

5.1051

5.8905

x =

0.3827

0.9239

0.3827

-0.3827

-0.9239

-0.3827

y =

0.9239

0.3827

-0.3827

-0.9239

-0.3827

0.3827

0.9239

>> fill(x,y,’r’) %fill creates a filled polygon based on x and y ; here points of x and y are 8 points on a circle circumference

>> axis square off, text(0,0,’STOP’,’Color’,[1 1 1],’FontSize’,80,’FontWeight’,’bold’,’HorizontalAlignment’,’center’)

>> a=[0.5 1 1.6 1.2 0.8 2.1]

a =

0.5000 1.0000 1.6000 1.2000 0.8000 2.1000

>> pie(a,a==max(a)) %create pie chart and take out the maximum

>> x=-2*pi:pi/10:2*pi, y=sin(x), z=3*cos(x)

x =

Columns 1 through 12

-6.2832 -5.9690 -5.6549 -5.3407 -5.0265 -4.7124 -4.3982 -4.0841 -3.7699 -3.4558 -3.1416 -2.8274

Columns 13 through 24

-2.5133 -2.1991 -1.8850 -1.5708 -1.2566 -0.9425 -0.6283 -0.3142 0 0.3142 0.6283 0.9425

Columns 25 through 36

1.2566 1.5708 1.8850 2.1991 2.5133 2.8274 3.1416 3.4558 3.7699 4.0841 4.3982 4.7124

Columns 37 through 41

5.0265 5.3407 5.6549 5.9690 6.2832

y =

Columns 1 through 12

0.0000 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 -0.0000 -0.3090

Columns 13 through 24

-0.5878 -0.8090 -0.9511 -1.0000 -0.9511 -0.8090 -0.5878 -0.3090 0 0.3090 0.5878 0.8090

Columns 25 through 36

0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000 -0.3090 -0.5878 -0.8090 -0.9511 -1.0000

Columns 37 through 41

-0.9511 -0.8090 -0.5878 -0.3090 -0.0000

z =

Columns 1 through 12

3.0000 2.8532 2.4271 1.7634 0.9271 -0.0000 -0.9271 -1.7634 -2.4271 -2.8532 -3.0000 -2.8532

Columns 13 through 24

-2.4271 -1.7634 -0.9271 0.0000 0.9271 1.7634 2.4271 2.8532 3.0000 2.8532 2.4271 1.7634

Columns 25 through 36

0.9271 0.0000 -0.9271 -1.7634 -2.4271 -2.8532 -3.0000 -2.8532 -2.4271 -1.7634 -0.9271 -0.0000

Columns 37 through 41

0.9271 1.7634 2.4271 2.8532 3.0000

>> subplot(2,1,1), plot(x,y,x,z)

>> subplot(2,1,2), plotyy(x,y,x,z) %plotyy causes the second plot to have different y-axis

>> x=-2.9:0.2:2.9, y=exp(-x.*x)

x =

Columns 1 through 12

-2.9000 -2.7000 -2.5000 -2.3000 -2.1000 -1.9000 -1.7000 -1.5000 -1.3000 -1.1000 -0.9000 -0.7000

Columns 13 through 24

-0.5000 -0.3000 -0.1000 0.1000 0.3000 0.5000 0.7000 0.9000 1.1000 1.3000 1.5000 1.7000

Columns 25 through 30

1.9000 2.1000 2.3000 2.5000 2.7000 2.9000

y =

Columns 1 through 12

0.0002 0.0007 0.0019 0.0050 0.0122 0.0271 0.0556 0.1054 0.1845 0.2982 0.4449 0.6126

Columns 13 through 24

0.7788 0.9139 0.9900 0.9900 0.9139 0.7788 0.6126 0.4449 0.2982 0.1845 0.1054 0.0556

Columns 25 through 30

0.0271 0.0122 0.0050 0.0019 0.0007 0.0002

>> subplot(2,2,1), bar(x,y) %bar plot

>> subplot(2,2,2), bar3(x,y) %in 3 D

>> subplot(2,2,3), stairs(x,y)

>> subplot(2,2,4), barh(x,y)

>> x=-2.9:0.2:2.9, y=randn(5000,1)

x =

Columns 1 through 11

-2.9000 -2.7000 -2.5000 -2.3000 -2.1000 -1.9000 -1.7000 -1.5000 -1.3000 -1.1000 -0.9000

Columns 23 through 30

1.5000 1.7000 1.9000 2.1000 2.3000 2.5000 2.7000 2.9000

y =

0.5377

0.3188

0.6762

0.5466

0.1106

…~~

-0.3271

1.3529

>> hist(y,x) %histogram plot

>> z=randn(30,1)

z =

-0.4277

-0.5794

-0.8450

-0.8051

>> stem(z,’–‘) %plot discrete sequence data; z may be matrix (then it would plot each column as stem plot)

>> set(gca,’YGrid’,’on’)

>> x=linspace(0,2,21), y=erf(x),e=rand(size(x))/10 %erf calculates error function (2*integration(exp(-t^2))~0-x~/sqrt(pi))

x =

Columns 1 through 12

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000

Columns 13 through 21

1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000

y =

Columns 1 through 12

0 0.1125 0.2227 0.3286 0.4284 0.5205 0.6039 0.6778 0.7421 0.7969 0.8427 0.8802

Columns 13 through 21

0.9103 0.9340 0.9523 0.9661 0.9763 0.9838 0.9891 0.9928 0.9953

e =

Columns 1 through 12

0.0345 0.0329 0.0927 0.0756 0.0288 0.0606 0.0766 0.0846 0.0902 0.0596 0.0069 0.0218

Columns 13 through 21

0.0869 0.0414 0.0661 0.0783 0.0248 0.0554 0.0230 0.0007 0.0767

>> errorbar(x,y,e) %plot error bar; here it plots x-vs-y and e(i) determines the length of the bar at each y(i)

>> t=linspace(0,2*pi),r=sin(2*t).*cos(2*t)

t =

Columns 1 through 12

0 0.0635 0.1269 0.1904 0.2539 0.3173 0.3808 0.4443 0.5077 0.5712 0.6347 0.6981

Columns 97 through 100

6.0928 6.1563 6.2197 6.2832

r =

Columns 1 through 12

0 0.1256 0.2431 0.3450 0.4249 0.4775 0.4994 0.4894 0.4480 0.3779 0.2835 0.1710

Columns 97 through 100

-0.3450 -0.2431 -0.1256 -0.0000

>> subplot(2,2,1), polar(t,r), title(‘Polar Plot’) %polar plot

>> z=eig(rand(20))

z =

9.6775 + 0.0000i

1.0852 + 0.5189i

-0.0587 + 0.0000i

0.0862 + 0.0000i

>> subplot(2,2,2), compass(z), title(‘Complex Number: Compass Plot’) %complex number plot; magnitude & angle, centered at origin

>> subplot(2,2,3), feather(z), title(‘Complex Number: Feather Plot’) %here, it plots each number individually at equal distances; at x=0,1,2. … 20

>> v=randn(1000,1)*pi

v =

3.3515

…~~

-0.2988

-3.1259

>> subplot(2,2,4), rose(v), title(‘Polar Histogram’) %polar histogram plot

>> x=rand(40,1), y=randn(40,1), area=20+(1:40)

x =

0.4944

0.3755

y =

-0.5401

0.0359

area =

Columns 1 through 20

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

Columns 21 through 40

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

>> scatter(x,y,area), box on %scatter plot x-vs-y and area determines the area of each circle

11.11 Easy Plotting

>> fplot(@humps,[-.5 3]), xlabel(‘x’), ylabel(‘hump(x)’), title(‘Humps function plot’) %fplot just takes function name and limits; here the function is humps function

>> fstr=’sin(x)/(x)’

fstr =

sin(x)/(x)

>> ezplot(fstr,[-15,15]) %ezplot takes the name of the function (which was defined earlier symbolically) and x range

>> istr='(x-2)^2/(2^2)+(y+1)^2/(3^2)-1′

istr =

(x-2)^2/(2^2)+(y+1)^2/(3^2)-1

>> ezplot(istr,[-2 6 -5 3]), axis square, grid %here y-axis limits is given too

11.12 Text Formatting

>> axis([0 1 0 0.5])

>> xlabel({‘first line’,’second line’}) %that’s how we can make multiple lines

>> text(0.2,0.1,’\itE=M\cdotC^{\rm2}’) %places text at (0.2,0.1); \it makes it italic

>> text(0.2,0.2,’\fontsize{16} \nabla \times H = J + \partialD/\partialt’) %placing something between {} is TEX

>> text(0.2,0.3,’\fontname{courier} \fontsize{16} \bf x_{\alpha} + y^{2\pi}’)

>> f=’f(t)=A_o+\fontsize{30}_\Sigma\fontsize{10}’

f =

f(t)=A_o+\fontsize{30}_\Sigma\fontsize{10}

>> text(0.2,0.4,[f ‘[A_ncos(n\omega_ot)+B_nsin(n\omega_ot)]’])

12. 3D Graphics

12.1 Line Plots

>> t=linspace(0,10*pi)

t =

Columns 1 through 11

0 0.3173 0.6347 0.9520 1.2693 1.5867 1.9040 2.2213 2.5387 2.8560 3.1733

Column 100

31.4159

>> plot3(sin(t), cos(t), t) %that’s a 3D plot for sine and cosine ; helical

>> xlabel(‘sin(t)’), ylabel(‘cos(t)’), zlabel(‘t’), grid on, title(‘sine and cosine in 3D’), text(0,0,0,’Origin’), v=axis %like before

v =

-1 1 -1 1 0 40

>> x=linspace(0,3*pi)

x =

Columns 1 through 12

0 0.0952 0.1904 0.2856 0.3808 0.4760 0.5712 0.6664 0.7616 0.8568 0.9520 1.0472

Columns 97 through 100

9.1392 9.2344 9.3296 9.4248

>> z1=sin(x),z2=sin(2*x),z3=sin(3*x), y1=zeros(size(x)), y3=ones(size(x)), y2=y3/2

z1 =

Columns 1 through 11

0 0.0951 0.1893 0.2817 0.3717 0.4582 0.5406 0.6182 0.6901 0.7557 0.8146

Column 100

0.0000

z2 =

Columns 1 through 11

0 0.1893 0.3717 0.5406 0.6901 0.8146 0.9096 0.9718 0.9989 0.9898 0.9450

Column 100

-0.0000

z3 =

Columns 1 through 11

0 0.2817 0.5406 0.7557 0.9096 0.9898 0.9898 0.9096 0.7557 0.5406 0.2817

Column 100

0.0000

y1 =

Columns 1 through 19

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Columns 96 through 100

0 0 0 0 0

y3 =

Columns 1 through 19

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

Columns 96 through 100

1 1 1 1 1

y2 =

Columns 1 through 11

0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000

Column 100

0.5000

>> plot3(x,y1,z1,x,y2,z2,x,y3,z3), grid on, xlabel(‘X’), ylabel(‘Y’),zlabel(‘Z’), title(‘sin(x), sin(2x), sin(3x)’)

>> plot3(x,z1,y1,x,z2,y2,x,z3,y3), grid on, xlabel(‘X’), ylabel(‘Y’),zlabel(‘Z’), title(‘sin(x), sin(2x), sin(3x)’)

12.2 Scalar Functions of 2 Variables

>> x=-3:3, y=1:5

x =

-3 -2 -1 0 1 2 3

y =

1 2 3 4 5

>> [X,Y]=meshgrid(x,y) %expands x and y as a mesh as seen below; X expand 5 rows and Y expand 6 columns (length y,x)

X =

-3 -2 -1 0 1 2 3

Y =

1 1 1 1 1 1 1

2 2 2 2 2 2 2

3 3 3 3 3 3 3

4 4 4 4 4 4 4

5 5 5 5 5 5 5

>> Z=(X+Y).^2

Z =

4 1 0 1 4 9 16

1 0 1 4 9 16 25

0 1 4 9 16 25 36

1 4 9 16 25 36 49

4 9 16 25 36 49 64

12.3 Mesh Plots

>> [X,Y,Z]=peaks(30) %generate 30 X 30 matrix of translating & scaling Gaussian distribution

X =

Columns 1 through 12

-3.0000 -2.7931 -2.5862 -2.3793 -2.1724 -1.9655 -1.7586 -1.5517 -1.3448 -1.1379 -0.9310 -0.7241

Columns 25 through 30

1.9655 2.1724 2.3793 2.5862 2.7931 3.0000

Y =

Columns 1 through 12

-3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -3.0000

Columns 25 through 30

3.0000 3.0000 3.0000 3.0000 3.0000 3.0000

Z =

Columns 1 through 12

0.0001 0.0002 0.0005 0.0011 0.0021 0.0036 0.0051 0.0048 -0.0001 -0.0138 -0.0403 -0.0810

Columns 25 through 30

0.0065 0.0028 0.0011 0.0004 0.0001 0.0000

>> mesh(X,Y,Z) %plot X, Y, Z in 3D

>> [X,Y,Z]=sphere(12) %generate 13 X 13 matrix of sphere coordinates

X =

Columns 1 through 12

0 0 0 0 0 0 0 0 0 0 0 0

-0.2588 -0.2241 -0.1294 0.0000 0.1294 0.2241 0.2588 0.2241 0.1294 0.0000 -0.1294 -0.2241

0 0 0 0 0 0 0 0 0 0 0 0

Column 13

-0.2588

Y =

Columns 1 through 12

0 0 0 0 0 0 0 0 0 0 0 0

0 -0.1294 -0.2241 -0.2588 -0.2241 -0.1294 0 0.1294 0.2241 0.2588 0.2241 0.1294

0 0 0 0 0 0 0 0 0 0 0 0

Column 13

Z =

Columns 1 through 12

-1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -1.0000

Column 13

1.0000

>> subplot(1,2,1), mesh(X,Y,Z), title(‘Opaque’), axis square off, hidden on %hidden on makes it opaque

>> subplot(1,2,2), mesh(X,Y,Z), title(‘Transparent’), axis square off, hidden off %hidden off makes it transparent

>> [X,Y,Z]=peaks(30);

>> meshc(X,Y,Z), title(‘Mesh plot with contours’) %meshc adds contour plots to the figure

>> meshz(X,Y,Z), title(‘Mesh plot with zero plane’)

>> waterfall(X,Y,Z), title(‘Waterfall plot’)

12.4 Surface Plots

>> [X,Y,Z]=peaks(30);

>> surf(X,Y,Z)

>> surf(X,Y,Z), shading flat

>> surf(X,Y,Z), shading interp %shading also have parametes pcolor and fill

>> x=X(1,:),y=Y(:,1),i=find(y>.8&y>1.2),j=find(x>-.6&x<.5),Z(i,j)=nan; %remove this region (0.88<y<1.2 & -0.6<x<0.5), thus not plot

>> surf(X,Y,Z)

>> [X,Y,Z]=peaks(30);

>> surfc(X,Y,Z) %add contour plot

>> surfl(X,Y,Z), shading interp, colormap pink

>> surfnorm(X,Y,Z) %add normals to the surface

>> [Nx,Ny,Nz]=surfnorm(X,Y,Z)

Nx =

Columns 1 through 11

-0.0002 -0.0010 -0.0021 -0.0040 -0.0062 -0.0071 -0.0028 0.0125 0.0449 0.0964 0.1582

Columns 12 through 22

0.2085 0.2226 0.1910 0.1219 0.0324 -0.0601 -0.1389 -0.1894 -0.2018 -0.1780 -0.1330

-0.2073 -0.1751 -0.1159 -0.0403 0.0405 0.1157 0.1744 0.2062 0.2052 0.1753 0.1303

Columns 23 through 30

-0.0857 -0.0484 -0.0241 -0.0107 -0.0042 -0.0014 -0.0004 0.0001

0.0855 0.0504 0.0268 0.0130 0.0058 0.0023 0.0009 0.0000

Ny =

Columns 1 through 11

-0.0002 -0.0006 -0.0015 -0.0034 -0.0066 -0.0108 -0.0137 -0.0094 0.0128 0.0666 0.1621

0.0001 0.0004 0.0013 0.0037 0.0095 0.0223 0.0482 0.0952 0.1707 0.2762 0.4005

Columns 12 through 22

0.2938 0.4357 0.5568 0.6381 0.6739 0.6644 0.6115 0.5197 0.4006 0.2753 0.1673

0.5226 0.6230 0.6915 0.7256 0.7256 0.6915 0.6231 0.5228 0.4008 0.2765 0.1711

Columns 23 through 30

0.0903 0.0435 0.0188 0.0073 0.0025 0.0008 0.0002 0.0000

0.0955 0.0485 0.0225 0.0096 0.0037 0.0013 0.0004 0.0001

Nz =

Columns 1 through 11

0.9983 0.9795 0.9205 0.7874 0.6135 0.4674 0.3667 0.2905 0.2214 0.1667 0.1343

1.0000 1.0000 1.0000 1.0000 0.9999 0.9994 0.9976 0.9918 0.9767 0.9449 0.8928

Columns 12 through 22

0.9328 0.8721 0.8084 0.7602 0.7381 0.7449 0.7790 0.8331 0.8937 0.9447 0.9769

0.8270 0.7624 0.7130 0.6869 0.6869 0.7131 0.7625 0.8271 0.8929 0.9449 0.9766

Columns 23 through 30

0.5941 0.8141 0.9486 0.9907 0.9988 0.9999 1.0000 1.0000

0.9917 0.9976 0.9994 0.9999 1.0000 1.0000 1.0000 1.0000

12.5 Mesh and Surface Plots for Irregular Data

>> x=rand(1,50), y=rand(1,50),z=peaks(x,y*pi)

x =

Columns 1 through 12

0.8147 0.9058 0.1270 0.9134 0.6324 0.0975 0.2785 0.5469 0.9575 0.9649 0.1576 0.9706

Columns 49 through 50

0.7094 0.7547

y =

Columns 1 through 12

0.2760 0.6797 0.6551 0.1626 0.1190 0.4984 0.9597 0.3404 0.5853 0.2238 0.7513 0.2551

Columns 49 through 50

0.0119 0.3371

z =

Columns 1 through 12

2.1033 2.0711 5.2532 2.0484 0.8004 8.0099 0.2598 3.4291 2.9515 2.1047 2.7197 2.1515

Columns 49 through 50

1.3328 2.9662

>> t=delaunay(x,y) %extending delaunary to 3D

t =

47 33 17

…

9 13 12

>> trimesh(t,x,y,z), hidden off, title(‘Triangular Mesh Plot’) %extending the triangular mesh plot to 3D

>> trisurf(t,x,y,z), title(‘Triangualr Surface Plot’)

>> voronoi(x,y,t) %extending to 3D

12.6 Changing Viewpoints

>> x=-7.5:.5:7.5, y=x

x =

Columns 1 through 12

-7.5000 -7.0000 -6.5000 -6.0000 -5.5000 -5.0000 -4.5000 -4.0000 -3.5000 -3.0000 -2.5000 -2.0000

Columns 25 through 31

4.5000 5.0000 5.5000 6.0000 6.5000 7.0000 7.5000

y =

Columns 1 through 12

-7.5000 -7.0000 -6.5000 -6.0000 -5.5000 -5.0000 -4.5000 -4.0000 -3.5000 -3.0000 -2.5000 -2.0000

Columns 25 through 31

4.5000 5.0000 5.5000 6.0000 6.5000 7.0000 7.5000

>> [X,Y]=meshgrid(x,y)

X =

Columns 1 through 11

-7.5000 -7.0000 -6.5000 -6.0000 -5.5000 -5.0000 -4.5000 -4.0000 -3.5000 -3.0000 -2.5000

Columns 12 through 22

-2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000

Columns 23 through 31

3.5000 4.0000 4.5000 5.0000 5.5000 6.0000 6.5000 7.0000 7.5000

Y =

Columns 1 through 11

Columns 12 through 22

Columns 23 through 31

>> R=sqrt(X.^2+Y.^2)+eps

R =

Columns 1 through 11

Columns 12 through 22

Columns 23 through 31

>> Z=sin(R)./R

Z =

Columns 1 through 11

Columns 12 through 22

Columns 23 through 31

>> subplot(2,2,1), surf(X,Y,Z), view(-37.5,30), title(‘Azimuth angle = -37.5 and elevation = 30’)

>> subplot(2,2,2), surf(X,Y,Z), view(-37.5+90,30), title(‘Azimuth angle = 52.5 and elevation = 30’)

>> subplot(2,2,3), surf(X,Y,Z), view(-37.5,60), title(‘Azimuth angle = -37.5 and elevation = 60’)

>> subplot(2,2,4), surf(X,Y,Z), view(0,90), title(‘Azimuth angle = 0 and elevation = 90’)

>> help view

view 3-D graph viewpoint specification.

view(AZ,EL) and view([AZ,EL]) set the angle of the view from which an

observer sees the current 3-D plot. AZ is the azimuth or horizontal

rotation and EL is the vertical elevation (both in degrees). Azimuth

revolves about the z-axis, with positive values indicating counter-

clockwise rotation of the viewpoint. Positive values of elevation

correspond to moving above the object; negative values move below.

view([X Y Z]) sets the view angle in Cartesian coordinates. The

magnitude of vector X,Y,Z is ignored.

Here are some examples:

AZ = -37.5, EL = 30 is the default 3-D view.

AZ = 0, EL = 90 is directly overhead and the default 2-D view.

AZ = EL = 0 looks directly up the first column of the matrix.

AZ = 180 is behind the matrix.

view(2) sets the default 2-D view, AZ = 0, EL = 90.

view(3) sets the default 3-D view, AZ = -37.5, EL = 30.

[AZ,EL] = view returns the current azimuth and elevation.

T = view returns the current general 4-by-4 transformation matrix.

view(AX,…) uses axes AX instead of the current axes.

12.7 Camera Control

12.8 Contour Plots

>> [X,Y,Z]=peaks();

>> contour(X,Y,Z,20)

>> contour3(X,Y,Z,20)

>> pcolor(X,Y,Z), shading interp

>> contourf(X,Y,Z,12) %filled contour lines

>> C=contour(X,Y,Z,12), clabel(C) %adding labels to the contour lines

>> [C,h]=contour(X,Y,Z,8), clabel(C,h)

12.9 Specialized 3D Plots

>> Z=peaks;

>> ribbon(Z)

>> Z=peaks;

>> ribbon(Z)

>> [X,Y,Z]=peaks(16);

>> [DX,DY]=gradient(Z,.5,.5); %find the gradient with space between points as .5 in x direction and .5 in y direction

>> quiver(X,Y,DX,DY) %plot direction field (gradients)

>> contour(X,Y,Z,10), hold on, quiver(X,Y,DX,DY)

>> [X,Y,Z]=peaks(20);

>> [Nx,Ny,Nz]=surfnorm(X,Y,Z);

>> surf(X,Y,Z), hold on, quiver3(X,Y,Z,Nx,Ny,Nz)

>> fill3(rand(3,5),rand(3,5),rand(3,5),rand(3,5)), grid on

>> Z=rand(5);

>> stem3(Z,’ro’,’filled’), grid on

12.10 Volume Visualization

>> x=linspace(-3,3,13); y=1:20; z=-5:5; [X,Y,Z]=meshgrid(x,y,z); size(X)

ans =

20 13 11

>> V=sqrt(X.^2+cos(Y).^2+Z.^2);

>> slice(X,Y,Z,V,[0 3],[5 15],[-3 5])

>> [xs,ys]=meshgrid(x,y);

>> zs=sin(-xs+ys/2);

>> slice(X,Y,Z,V,xs,ys,zs)

>> slice(X,Y,Z,V,[0 3],[5 15],[-3 5]), h=contourslice(X,Y,Z,V,3,[5 15],[]), set(h,’EdgeColor’,’k’,’Linewidth’,1.5)

>> [X,Y,Z,V]=flow(13); fv=isosurface(X,Y,Z,V,-2);

>> subplot(1,2,1), p=patch(fv); set(p,’FaceColor’,[.5 .5 .5],’EdgeColor’,’Black’), view(3), axis equal tight, grid on

>> subplot(1,2,2), p=patch(shrinkfaces(fv,.3)); set(p,’FaceColor’,[.5 .5 .5],’EdgeColor’,’Black’), view(3), axis equal tight, grid on

>> [X,Y,Z,V]=flow; fv=isosurface(X,Y,Z,V,-2);

>> subplot(2,2,1), p=patch(fv); Np=size(get(p,’Faces’),1) ,set(p,’FaceColor’,[.5 .5 .5],’EdgeColor’,’Black’), view(3), axis equal tight, grid on

Np =

3568

>> subplot(2,2,2),[Xr,Yr,Zr,Vr]=reducevolume(X,Y,Z,V,[3 2 2]); fvr=isosurface(Xr,Yr,Zr,Vr,-2); p=patch(fvr); Np=size(get(p,’Faces’),1) ,set(p,’FaceColor’,[.5 .5 .5],’EdgeColor’,’Black’), view(3), axis equal tight, grid on

Np =

656

>> subplot(2,2,3), p=patch(fv); set(p,’FaceColor’,[.5 .5 .5],’EdgeColor’,’Black’), view(3), axis equal tight, grid on, reducepatch(p,.15),Np=size(get(p,’Faces’),1),

Np =

534

>> subplot(2,2,4), p=patch(fvr); set(p,’FaceColor’,[.5 .5 .5],’EdgeColor’,’Black’), view(3), axis equal tight, grid on, reducepatch(p,.15),Np=size(get(p,’Faces’),1),

Np =

>> data=rand(10,10,10);

>> datas=smooth3(data,’box’,3);

>> subplot(1,2,1)

>> p=patch(isosurface(data,.5),’FaceColor’,’Blue’,’EdgeColor’,’none’); patch(isocaps(data,.5),’FaceColor’,’interp’,’EdgeColor’,’none’); isonormals(data,p), view(3); axis vis3d tight off; camlight; lighting phong

>> subplot(1,2,2)

>> p=patch(isosurface(datas,.5),’FaceColor’,’Blue’,’EdgeColor’,’none’); patch(isocaps(datas,.5),’FaceColor’,’interp’,’EdgeColor’,’none’); isonormals(datas,p), view(3); axis vis3d tight off; camlight; lighting phong

12.11 Easy Plotting

>> f=[‘3*(1-x).^2.*exp(-(x.^2)-(y+1).^2)-10*(x/5-x.^3-y.^5).*exp(-x.^2-y.^2)-1/3*exp(-(x+1).^2-y.^2)’]

f =

3*(1-x).^2.*exp(-(x.^2)-(y+1).^2)-10*(x/5-x.^3-y.^5).*exp(-x.^2-y.^2)-1/3*exp(-(x+1).^2-y.^2)