//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 =
1
2
3
4
5
6
7
>> 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 =
1
2
3
4
>> 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
1 1 1 1
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
1 1 1 1
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
1 1 1
1 1 1
ans =
0 0 0 0 0
0 0 0 0 0
ans =
1 1 1 1
1 1 1 1
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
1 1 1
1 1 1
ans =
0 0 0 0 0
0 0 0 0 0
ans =
1 1 1 1
1 1 1 1
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
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
3.1416 3.1416 3.1416 3.1416
3.1416 3.1416 3.1416 3.1416
ans =
3.1416 3.1416 3.1416 3.1416
3.1416 3.1416 3.1416 3.1416
3.1416 3.1416 3.1416 3.1416
ans =
3.1416 3.1416 3.1416 3.1416
3.1416 3.1416 3.1416 3.1416
3.1416 3.1416 3.1416 3.1416
ans =
3.1416 3.1416 3.1416 3.1416
3.1416 3.1416 3.1416 3.1416
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 =
2
4
8
>> A(:,[2]) %same as above
ans =
2
4
8
>> A(:,[2 2]) %this means 2nd column 2 times
ans =
2 2
4 4
8 8
>> B=A(:,[2 2 2 2]) %create B matrix by repeating 2nd 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 =
1
1
1
1
1
1
>> H(:)=A(:,2:3) %you guess what happens here 😉
H =
2
4
8
3
7
9
>> 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
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 =
2
>> D(5)
ans =
5
>> D(end)
ans =
12
>> 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 =
11
>> [r,c]=ind2sub(size(D),11) %the converse
r =
2
c =
4
>> 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
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 =
1
1
3
6
idx =
1
2
4
3
>> 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
1 2 3 4 5 6
1 2 3 4 5 6
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 =
3
3
2
3
c =
1
2
3
3
>> k=find(A>5) %returns single indices. Again, remember it searches column by column
k =
3
6
8
9
>> A(k) %get the elements addressed by k, i.e., those > 5
ans =
7
8
6
9
>> 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
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 =
7
8
6
9
>> B(k) %just to check
ans =
7
8
6
9
>> B(r,c)=0 %indexing like above and setting the resulted elements to 0
B =
1 2 3
0 0 0
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 =
4
>> [mx,idm]=min(v) %minimum value and its location
mx =
0.0975
idm =
6
>> 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 =
6
>> 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 =
13
>> 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 =
1
5
9
>> 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
3.1416 3.1416 3.1416 3.1416
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 =
2
c =
3
>> r=size(A,1) %here we ask only for number of rows, the argument 1 tells MATLAB so
r =
2
>> c=size(A,2) %here we ask only for number of columns, the argument 2 tells MATLAB so
c =
3
>> numel(A) %returns total number of elements in matrix A
ans =
6
>> length(A) %returns number of rows or columns, whichever is larger
ans =
3
>> 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 =
0
>> max(size(d)) %it would make sense to use length(d) instead, but d has a zero dimension!
ans =
3
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
0 0 0
0 0 0
0 0 0
A(:,:,2) =
0 0 0
0 0 0
0 0 0
0 0 0
>> A=zeros(2,3) %that’s a 2D array
A =
0 0 0
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
0 0 0
A(:,:,2) =
1 1 1
1 1 1
>> A(:,:,3)=4 %adding a third page
A(:,:,1) =
0 0 0
0 0 0
A(:,:,2) =
1 1 1
1 1 1
A(:,:,3) =
4 4 4
4 4 4
>> B=reshape(A,2,9) %create 2D array
B =
0 0 0 1 1 1 4 4 4
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
0 0 0 1 1 1 4 4 4
>> reshape(B,2,3,3) %turn into 3D array
ans(:,:,1) =
0 0 0
0 0 0
ans(:,:,2) =
1 1 1
1 1 1
ans(:,:,3) =
4 4 4
4 4 4
>> reshape(B,[2 3 3]) %same as above
ans(:,:,1) =
0 0 0
0 0 0
ans(:,:,2) =
1 1 1
1 1 1
ans(:,:,3) =
4 4 4
4 4 4
>> C=ones(2,3)
C =
1 1 1
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
1 1 1
ans(:,:,2) =
1 1 1
1 1 1
ans(:,:,3) =
1 1 1
1 1 1
>> a=zeros(2)
a =
0 0
0 0
>> b=ones(2)
b =
1 1
1 1
>> c=repmat(2,2,2)
c =
2 2
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
0 0
D(:,:,2) =
1 1
1 1
D(:,:,3) =
2 2
2 2
>> E=cat(4,a,b,c) %here 4D array from a, b, c
E(:,:,1,1) =
0 0
0 0
E(:,:,1,2) =
1 1
1 1
E(:,:,1,3) =
2 2
2 2
>> E(:,1,:,:) %see 1st column in all rows, all pages, all 4th dimension
ans(:,:,1,1) =
0
0
ans(:,:,1,2) =
1
1
ans(:,:,1,3) =
2
2
>> 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
0 0
F(:,:,2) =
1 1
1 1
F(:,:,3) =
2 2
2 2
>> size(F)
ans =
2 2 3
>> v(1,1,:)=1:6 %page-vector
v(:,:,1) =
1
v(:,:,2) =
2
v(:,:,3) =
3
v(:,:,4) =
4
v(:,:,5) =
5
v(:,:,6) =
6
>> squeeze(v) %turns into column vector
ans =
1
2
3
4
5
6
>> v(:) %same
ans =
1
2
3
4
5
6
>> G=cat(3,2+zeros(2,4),ones(2,4),zeros(2,4))
G(:,:,1) =
2 2 2 2
2 2 2 2
G(:,:,2) =
1 1 1 1
1 1 1 1
G(:,:,3) =
0 0 0 0
0 0 0 0
>> H=reshape(G,[3 2 4])
H(:,:,1) =
2 2
2 2
2 2
H(:,:,2) =
2 1
2 1
1 1
H(:,:,3) =
1 1
1 0
1 0
H(:,:,4) =
0 0
0 0
0 0
>> K=reshape(G,[4 3 2])
K(:,:,1) =
2 2 1
2 2 1
2 2 1
2 2 1
K(:,:,2) =
1 0 0
1 0 0
1 0 0
1 0 0
>> L=reshape(G,2,12)
L =
2 2 2 2 1 1 1 1 0 0 0 0
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 =
1
>> sub2ind(size(G),1,2,1) %find 1st row, 2nd column, 1st page element
ans =
3
>> sub2ind(size(G),1,2,3) %find 1st row, 2nd column, 3rd page element
ans =
19
>> [r,c,p]=ind2sub(size(G),19) %find where elemet 19 is located
r =
1
c =
2
p =
3
>> 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 =
18
>> [r,c,p]=size(M)
r =
2
c =
3
p =
3
>> r=size(M,1)
r =
2
>> c=size(M,2)
c =
3
>> p=size(M,3)
p =
3
>> v=size(M,4)
v =
1
>> nm=size(M,5)
nm =
1
>> nm=size(M,13)
nm =
1
>> nm=size(M,50)
nm =
1
>> ndims(M)
ans =
3
>> ndims(M(:,:,1))
ans =
2
>> length(size(M)) %equivalent to ndims(M)
ans =
3
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 1st and 3rd 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 3rd 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]
[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]
[ 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 2nd row 2nd column and put it in variable x
x =
5
>> y=B(2,2) %take the whole cell located in 2nd row 2nd 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 =
1
ans =
0
>> isa(y,’cell’) %check if y is a cell
ans =
1
>> isa(y{1},’double’) %check if content index 1 of y is a double
ans =
1
>> isa(y{1},’numeric’) %check if content index 1 of y is a numeric
ans =
1
>> isa(y{1},’cell’) ) %check if content index 1 of y is a cell
ans =
0
>> 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 =
1
>> B{:,2} %get all rows 2nd column
ans =
Hello World
ans =
5
>> 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
1 1 1
b =
0
0
c =
3
4
>> 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 =
2
cl =
5
>> 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]
[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
1 1 1
F{2} =
0
0
F{3} =
3
4
>> [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
1 1 1
s =
0
0
t =
3
4
>> [r,s,t]=deal(F{1},F{2},F{3}) %same as above
r =
1 1 1
1 1 1
s =
0
0
t =
3
4
>> [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 =
1
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]
[1×6 double]
[1×6 double]
2.6 Cell Arrays of Strings
>> s=’Hello World’
s =
Hello World
>> cs={‘Hi everyone’
s
‘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 =
1
>> sa=char(cs) %builds a char string array from z cell array
sa =
Hi everyone
Hello World
Welcome to planet Earth
>> ischar(sa)
ans =
1
>> iscell(sa)
ans =
0
>> 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 =
1
>> size(cst)
ans =
3 1
>> isequal(cs,cst)
ans =
1
>> isequal(cs,sa)
ans =
0
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 =
1
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 =
1
>> 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 =
y
>> circle(1).filled(2:end)
ans =
es
>> 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
1 1 1
b =
0
0
c =
3
4
>> 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 =
2
n =
5
>> 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 =
5
ans =
14
ans =
30
>> [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 =
1
>> isfield(circle,’width’)
ans =
0
>> class(square)
ans =
struct
>> isstruct(circle) %check if circle is a structure
ans =
1
>> isstruct(r1)
ans =
0
>> 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 =
5
>> r2=getfield(circle,{3},fnames{1})
r2 =
30
>> 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 =
30
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 y1 = x and y2 = 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 =
10
>> 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
N
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
N
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 =
5
>> 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 =
1
>> 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 =
1
>> 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 2nd 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
1
2
3
4
5
>> 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 =
53
>> EPS=EPS*2
EPS =
2.2204e-16
>> eps
ans =
2.2204e-16
7.3 If-Else-End Constructions
>> apples=10
apples =
10
>> cost=apples*25
cost =
250
>> if apples>5
cost=(1-20/100)*cost;
end
>> cost
cost =
200
>> EPS=1
EPS =
1
>> for num=1:100
EPS=EPS/2;
if(1+EPS)<=1
EPS=EPS*2
break
end
end
EPS =
2.2204e-16
>> num
num =
53
>> EPS=1
EPS =
1
>> 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
1 1
1 1
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
4 4
4 4
4 4
>> x=ones(4,2)
x =
1 1
1 1
1 1
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*************************
********************************M-FILE*************************
function v=myfun(varargin)
%varargin function is used when we need unlimited number of input arguments; 1st 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*************************
********************************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 =
3
y =
5
>> 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 =
45
>> Hz_min=Hzi(i) %find corresponding minimum Hz value
Hz_min =
3.3333e+03
>> x=linspace(0,2*pi,11)’
x =
0
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 =
20
j =
20
>> 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
1 2 3 4 5
1 2 3 4 5
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 =
2
>> DT(ans,:)
ans =
11 3 5
>> pointLocation(DT,[0,0;0,-0.5])
ans =
2
NaN
>> dsearchn([x’ y’],tri,[0 0.5; 1 2])
ans =
1
8
>> 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 =
0
>> [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 =
0
>> 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 =
0
>> 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 =
0
>> 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 =
0
>> options=optimset(‘Display’,’none’); %show nothing
>> [x,val]=fzero(@humps,[-2 0],options)
x =
-0.131618018099606
val =
0
>> 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’)
or
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 =
1
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 =
1
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 =
1
>> plot(x,z,’m’)
>> hold off
>> ishold
ans =
0
>> 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
1.1532
2.7661
m =
0.1842
0.9152
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.9239
0.3827
-0.3827
-0.9239
-0.9239
-0.3827
y =
0.9239
0.3827
-0.3827
-0.9239
-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
-3 -2 -1 0 1 2 3
-3 -2 -1 0 1 2 3
-3 -2 -1 0 1 2 3
-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.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
-0.2588
0
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
0
0
0
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 =
98
>> 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)
>> subplot(2,2,1), ezmesh(f)
>> subplot(2,2,2), ezsurf(f)
>> subplot(2,2,3), ezcontour(f)
>> subplot(2,2,4), ezcontourf(f)
Appendix
Common Tasks
Generate an array of real numbers
Generating an array of complex numbers
Generating x,y,x numbers
Plot, 2D
Plot, 3D