650 likes | 805 Views
math403 Introductory MATLAB. Zheng Chen Spring 2010. Environment of matlab. 3 windows: Command window History window Editor window. Executable file form. Files of form .m for functions and main programs All files to solve one problem should be in one directory
E N D
math403Introductory MATLAB Zheng Chen Spring 2010
Environment of matlab • 3 windows: • Command window • History window • Editor window
Executable file form • Files of form .m for functions and main programs • All files to solve one problem should be in one directory • Pwd, cd, ls, cd.., cp, mkdir, rmdir, • Help while
Matlab as a calculator • 2+3/4.0*5+(-2)^3+ log(3.4) • Sin(pi/4); • 1.2+3i • Format of numbers • Format long • Format short
variables • X=3-2^3; • Assign the value to x • Y=x+4; assign the value to y • Pi=3.1415…, eps=2.2204e-16 • I, j have the value sqrt(-1);
Built in functions • Trig • Log • Sqrt(x); • exp
vectors • Row vector: V=[1,2 3] • U=[-5, 6, 12] • V+/-U; • W=[7,8]; • Cd=[2*w, 3V]; • V(2); it is the same as V(1,2); • Column vector w=[7;8;9] • W(2) =w(2,1);
Colon notations • V=1:8; • V=3:7; • V=.32:0.1:-2 • Extracting • A=[1,2 ,3 ,4,5 ,6 ,7 ,8] • What is the b=A(4:7) and c=A(1:3:7)
transpose • Row vector to column vector or vice verse • General matrices with size m*n • W’ is of size n*m if W has size m*n
Scalar product of vectors u*v • Scalar product: • u=[a, b, c], v=[l, m,n]’ • In matlab, u*v=al+bm+cn; it is a number • sqrt(u*u’) ; and norm(u); • Recall what is the meaning of norm?
Dot product .* • u=[a, b, c], v=[l, m, n] • u.*v is a vector with the same size as u and v; • simply multiply the corresponding elements of two vectors • u.*v named as dot product
Dot Division ./ • ./ is defined to give element by element division • u=[a, b, c], v=[l, m, n]; • u./v=[a/l , b/m, c/n]; • Try 1./[2, 3, 4];
Dot Power of Arrays (.^) • u = [ 10, -11, 12]; • u.^2; • .* ; ./ ; .^ mean component wise operations
Functions applied to arrays • u = [ 10, 11, 12]; • sin(u), sqrt(u);
Matrices –Two-Dimensional Arrays • Input: a=[1 2 3; 5 6 7]; • or a=[1 2 3 5 6 7]; b=[1 2; 3 4; 5 6]; • size(a) = dimensions of matrix, that is, number of rows and number of columns • size(a) is itself one vector with size 1, 2.
Special Matrices • ones(2,3); zeros(2,3); eye(3); • d = [-3 4 2], D = diag(d); • given F = [0 1 8 7; 3 -2 -4 2; 4 2 1 1], find diag(F); • c=[0 1; 3 -2; 4 2]; x=[8;-4;1]; find g = [c x] (extra column) find h=[c; k] (stacked c and k) where k=[7 8] • Try spy(F), grid for fun
Extracting Bits of Matrices • J =[ 1 2 3 4 5 6 7 8 9 0 3 11]; Find J(2,4); J(:,2); J(3,:); J(1:2, 2:4);
Dot product of matrices (.*) • The dot product works as for vectors: corresponding elements are multiplied together
Matrix{vector products • An (m x n) matrix times an (n x k) matrix is a (m x k) matrix. • A = [5 7 9; 1 -3 -7]; x = [8; -4; 1]; find A*x • Practice: choose two matrices A and B, find A*B, (A*B)’ and B’*A’ • Practice: choose matrices A, B and C, find A*(B*C) and (A*B)*C
Sparse matrices • >> i = [1, 3, 5]; j = [2,3,4]; • >> v = [10 11 12]; • >> S = sparse (i,j,v); T=full(S); • save memory, only save diagonal elements • >> n = 5; • >> l = -(2:n+1)'; d = (1:n )'; u = ((n+1):-1:2)'; • >> B = spdiags([l d u],-1:1,n,n); • >> full(B) (%l,d,u are column vectors)
System of linear equations • Ax=b • x = A \ b; (“\backslash") • try inv(A) and inv(A)*b
Overdetermined system of linear equations • generally speaking, Ax=b has no solution when # of equations > # of unknowns • r = Ax-b; • Least square solution: make r minimum • normal equations: Cx = d; where C = A’A and d = A’b • matlab gives the solution in its routine library x = A\b
A =[0.0000 0.0000; 0.0001 0.0000; 0.0013 0.0000; 0.0030 0.0000; 0.0075 0.0001] • b = [5 50 500 1000 2000]; • X=A\b’
For Loop • for i=1:10 disp([i i^2 i^3]) end • for i=1:4 for j=1:4 disp([i j i^2+j^2]) end end
example • for i=1:4 for j=i:4 disp([i j i^2+j^2]) end end
example • sum=0.0; bigsum=0.0; • for i=i:1000 bigsum =bigsum+1/i; end for i=i:1000 sum =sum+1/i^2; end disp([bigsum sum])
Questions • What is the sum of 1+1/2+1/3+… until the term where 1/n<=10^(-3) ? • How many terms we have in the expansion of 1+1/2+1/3+… if the partial sum is just >=2 ?
While loop • % format long i = 1; % start the counter at 1 total = 0; % initialize the sum while(1/i >= 1e-3) % condition total = total + 1/i; % update the sum i=i+1; end disp([i 1/i total]) % display the number of loops, 1/i, and the current sum
While loop(condition changed) i = 1; % start the counter at 1 total = 0; % initialize the sum while(total <= 2) % condition changed total = total + 1/i % update the sum i= i+1 end % disp([I total]) disp([i-2 total-1/(i-1)]) % display the number of terms, and the final sum
Example to use while loop • Problem: what is the greatest value of n that can be used in the sum 1+(1/2)^2+ (½)^2 +… + (1/n)^2 and get a value of less than 2? S = 1; n = 1; while S+ (1/(n+1))^2 < 2 n = n+1; S = S + (1/n)^2; end • [n, S]
Logicals • true = 1, false = 0 • If x<1.0 elseif
if...then...else...end • if logical test 1 commands to be executed if test 1 is true elseif logical test 2 commands to be executed if test 2 is true but test 1 is false ... end
R = input('What is your name:','s') if R=='cuana' 'today your hw: read the notes' else 'bye' end
Problem: find out all triples of integers, at most 20, as the 3 sides of one right triangle, like 3, 4, 5;
Using for and if loop for i=1:20 for j=i:20 for k=j:20 if k^2==i^2+j^2 [i, j, k] end end end end
Functions defined by yourself • % Compute the area of a triangle whose • % sides have length a, b and c. • % Inputs: a,b,c: Lengths of sides • % Output: A: area of triangle • % Usage: Area = area(2,3,4); function [A] = area(a,b,c) s = (a+b+c)/2; A = sqrt(s*(s-a)*(s-b)*(s-c));
Problem: find the are of one triangle with given 3 numbers, you need to decide whether these three numbers can be the 3 sides of one triangle. for example, 3, 4 and 5 are 3 sides of one triangle, 3, 4 and 8 are not 3 sides of one triangle. given a<=b<=c , criteria for them to be 3 sides of a triangle: c<a+b & a> c-b.
Functions countinue • function [A] =triValue(a,b,c) s = (a+b+c)/2; A = sqrt(s*(s-a)*(s-b)*(s-c)); • save the file as triValue.m • Another function: two results are provided function [cir, A] =triValue(a,b,c) s = (a+b+c)/2; cir=s*2; A = sqrt(s*(s-a)*(s-b)*(s-c));
function area=tri_area(a, b, c) dd=min([a,b,c]); aa=max([a,b,c]); b=a+b+c-dd-aa; a=aa; c=dd; if c>(a-b) & b+c>a 'the given numbers are 3 sides of one triangle, and the area is: ' s=(a+b+c)/2; area=sqrt(s*(s-a)*(s-b)*(s-c)) else 'the given numbers are not 3 sides of one triangle' end
Another example of functions • Fibonnaci sequence is dened by n=3,4, … input: integrer output:
Method 1: File Fib1.m function f = Fib1(n) % Returns the nth number in the % Fibonacci sequence. F=zeros(1,n+1); F(2) = 1; for i = 3:n+1 F(i) = F(i-1) + F(i-2); end f = F(n);
function f = Fib2(n) % Returns the nth number in the Fibonacci sequence. if n==1 f = 0; elseif n==2 f = 1; else f1 = 0; f2 = 1; for i = 2:n-1 f = f1 + f2; f1=f2; f2 = f; end end
Recursive programming function f = Fib3(n) % Returns the nth number in the % Fibonacci sequence. if n==1 f = 0; elseif n==2 f = 1; else f = Fib3(n-1) + Fib3(n-2); end
% The final version uses matrix powers. The vector y has two components, % Returns the nth number in the Fibonacci sequence. function f = Fib4(n) A = [0 1;1 1]; y = A^n*[1;0]; f=y(1);
Which method you like most? • Which algorithm is most efficient? • method 3(recursive programming) no efficient
More Built-in Functions • x = pi*(-1:3), round(x) • Try fix(x); floor(x); ceil(x); sign(x), rem(x,3) • A = [1:3; 4:6; 7:9]; • s = sum(A), ss = sum(sum(A)); • x = [1.3 -2.4 0 2.3], max(x), max(max(A)); min(x); • y = rand, Y = rand(2,3)
Use rand % simulates "n" throws of a pair of dice % Input: n, the number of throws % Output: an n times 2 matrix, each row % referring to one throw. % Useage: T = dice(3) function [d] = dice(n) d = floor(1 + 6*rand(n,2)); %% end of dice
Function ‘find’ • A = [ -2 3 4 4; 0 5 -1 6; 6 8 0 1]; • k = find(A==0); To interpret the result we have to recognize That “find" first reshapes A into a column vector; this is equivalent to numbering the elements of A by columns as in 1 4 7 10 2 5 8 11 3 6 9 12 Now try n = find(A <= 0); and A(n)
A=[1 4 7 10 2 5 8 11 3 6 9 12]; n = find(A <= 0); Try A(n); Try m = find( A' == 0);
plot • linspace (a,b,n) producing n+1 points; • N = 10; h = 1/N; x = 0:h:1; doing the same. • X=linspace(0,1,110); • y = sin(3*pi*x); • Plot(x,y)