穨matlab教學範例ccc.doc

Size: px

Start display at page:

Download "穨matlab教學範例ccc.doc"

潮牛
5 years ago
Views:

1 MATLAB ( Math Dept. of University of New Hampshire) ( ) 1. Tutorial on vectors 2. Tutorial on matrices 3. Tutorial on vector operations 4. Tutorial on loops 5. Tutorial on plots 6. Tutorial on executable files 7. Tutorial on subroutines 8. Tutorial on if statements

2 >> 0:2: >> ans' >> v = [0:2:8] v = >> v v = >> v; >> v' >> v(1:3) 2

3 0 2 4 >> v(1:2:4) 0 4 >> v(1:2:4)' 0 4 >> v(1:3)-v(2:4) >> A = [ 1 2 3; 3 4 5; 6 7 8] A = >> B = [ [1 2 3]' [2 4 7]' [3 5 8]'] B = >> whos Name Size Elements Bytes Density Complex A 3 by Full No B 3 by Full No ans 1 by Full No v 1 by Full No 3

4 Grand total is 26 elements using 208 bytes >> v = [0:2:8] v = >> A*v(1:3) ??? Error using ==> * Inner matrix dimensions must agree. >> A*v(1:3)' >> A(1:2,3:4)??? Index exceeds matrix dimensions. >> A(1:2,2:3) >> A(1:2,2:3)' >> inv(a) Warning: Matrix is close to singular or badly scaled. 1.0e+15 * Results may be inaccurate. RCOND = e

5 >> inv(a)??? Undefined function or variable a. >> eig(a) >> [v,e] = eig(a) v = e = >> diag(e) >> v = [1 3 5]' v = 1 3 5

6 5 >> x = A\v Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = e-18 x = 1.0e+15 * >> x = B\v x = >> B*x >> x1 = v'/b x1 = >> x1*b >> clear >> whos >> v = [1 2 3]' v = 1 6

7 2 3 >> b = [2 4 6]' b = >> v+b >> v-b >> v*b Error using ==> * Inner matrix dimensions must agree. >> v*b' >> v'*b >> v.*b 2 7

8 8 18 >> v./b >> sin(v) >> log(v) >> x = [0:0.1:100] x = Columns 1 through (stuff deleted) Columns 995 through >> y = sin(x).*x./(1+cos(x)); >> plot(x,y) >> plot(x,y,'rx') >> help plot 8

9 PLOT Plot vectors or matrices. PLOT(X,Y) plots vector X versus vector Y. If X or Y is a matrix, then the vector is plotted versus the rows or columns of the matrix, whichever line up. PLOT(Y) plots the columns of Y versus their index. If Y is complex, PLOT(Y) is equivalent to PLOT(real(Y),imag(Y)). In all other uses of PLOT, the imaginary part is ignored. Various line types, plot symbols and colors may be obtained with PLOT(X,Y,S) where S is a 1, 2 or 3 character string made from the following characters: y yellow. point m magenta o circle c cyan x x-mark r red + plus g green - solid b blue * star w white : dotted k black -. dashdot -- dashed For example, PLOT(X,Y,'c+') plots a cyan plus at each data point. PLOT(X1,Y1,S1,X2,Y2,S2,X3,Y3,S3,...) combines the plots defined by the (X,Y,S) triples, where the X's and Y's are vectors or matrices and the S's are strings. For example, PLOT(X,Y,'y-',X,Y,'go') plots the data twice, with a solid yellow line interpolating green circles at the data points. The PLOT command, if no color is specified, makes automatic use of the colors specified by the axes ColorOrder property. The default ColorOrder is listed in the table above for color systems where the default is yellow for one line, and for multiple lines, to cycle through the first six colors in the table. For monochrome systems, PLOT cycles over the axes LineStyleOrder property. PLOT returns a column vector of handles to LINE objects, one handle per line. The X,Y pairs, or X,Y,S triples, can be followed by 9

10 parameter/value pairs to specify additional properties of the lines. See also SEMILOGX, SEMILOGY, LOGLOG, GRID, CLF, CLC, TITLE, XLABEL, YLABEL, AXIS, AXES, HOLD, and SUBPLOT. >> plot(x,y,'y',x,y,'go') >> plot(x,y,'y',x,y,'go',x,exp(x+1),'m--') >> whos Name Size Elements Bytes Density Complex ans 3 by Full No b 3 by Full No v 3 by Full No x 1 by Full No y 1 by Full No Grand total is 2011 elements using bytes >> coef = zeros(1,1001); >> coef(1) = y(1); >> y = (y(2:1001)-y(1:1000))./(x(2:1001)-x(1:1000)); >> whos Name Size Elements Bytes Density Complex ans 3 by Full No b 3 by Full No coef 1 by Full No v 3 by Full No x 1 by Full No y 1 by Full No Grand total is 3011 elements using bytes >> coef(2) = y(1); >> y(1) >> y = (y(2:1000)-y(1:999))./(x(3:1001)-x(1:999)); >> coef(3) = y(1); 10

11 >> for j=1:4, j j = 1 j = 2 j = 3 j = >> 4 >> v = [1:3:10] v = >> for j=1:4, v(j) = j; 11

12 >> v v = >> A = [ [1 2 3]' [3 2 1]' [2 1 3]'] A = >> B = A; >> for j=2:3, A = A = A(j,:) = A(j,:) - A(j-1,:) >> for j=2:3, B = B = for i=j:3, B(i,:) = B(i,:) - B(j-1,:)*B(i,j-1)/B(j-1,j-1)

13 B = >> h = 0.1; >> x = [0:h:2]; >> y = 0*x; >> y(1) = 1; >> size(x) 1 21 >> for i=2:21, y(i) = y(i-1) + h*(x(i-1)^2 - y(i-1)^2); >> plot(x,y) >> plot(x,y,'go') >> plot(x,y,'go',x,y) >> h = 0.001; >> x = [0:h:2]; >> y = 0*x; >> y(1) = 1; >> i = 1; >> size(x) >> max(size(x)) 2001 >> while(i<max(size(x))) y(i+1) = y(i) + h*(x(i)-abs(y(i))); 13

14 i = i + 1; >> plot(x,y,'go') >> plot(x,y) >> h = 1/16; >> x = 0:h:1; >> y = 0*x; >> size(y) 1 17 >> max(size(y)) 17 >> y(1) = 1; >> for i=2:max(size(y)), y(i) = y(i-1) + h/y(i-1); >> true = sqrt(2*x+1); >> plot(x,y,'go',x,true) >> plot(x,abs(true-y),'mx') >> subplot(1,2,1); >> plot(x,y,'go',x,true) >> subplot(1,2,2); >> plot(x,abs(true-y),'mx') 14

15 Figure 1. The two plots from the first approximation >> clf >> h = h/2; >> x1 = 0:h:1; >> y1 = 0*x1; >> y1(1) = 1; >> for i=2:max(size(y1)), y1(i) = y1(i-1) + h/y1(i-1); >> true1 = sqrt(2*x1+1); >> plot(x,y1,'go',x,true1)??? Error using ==> plot Vectors must be the same lengths. >> plot(x1,y1,'go',x1,true1) >> plot(x1,abs(true1-y1),'mx') >> subplot(1,2,1); >> plot(x,abs(true-y),'mx') >> subplot(1,2,2); >> plot(x1,abs(true1-y1),'mx') >> title('errors for h=1/32') 15

16 >> xlabel('x'); >> ylabel(' Error '); >> subplot(1,2,1); >> xlabel('x'); >> ylabel(' Error '); >> title('errors for h=1/16') Figure 2. The errors for the two approximations file: simpleeuler.m This matlab file will find the approximation to 16

17 dy/dx = 1/y y(0) = starty To run this file you will first need to specify the step the following: h : the step size starty : the initial value The routine will generate three vectors. The first vector is x which is the grid points starting at x0=0 and have a step size h. The second vector is an approximation to the specified D.E. The third vector is the true solution to the D.E. If you haven't guessed, you cna use the percent sign to add comments. x = [0:h:1]; y = 0*x; y(1) = starty; for i=2:max(size(y)), y(i) = y(i-1) + h/y(i-1); true = sqrt(2*x+1); >> simpleeuler??? Undefined function or variable h. Error in ==> /home/black/math/mat/examples/simpleeuler.m On line 27 ==> x = [0:h:1]; >> h = 1/16; 17

18 >> starty = 1; >> starty = 1; >> simpleeuler >> whos Name Size Elements Bytes Density Complex h 1 by Full No i 1 by Full No starty 1 by Full No true 1 by Full No x 1 by Full No y 1 by Full No Grand total is 54 elements using 432 bytes >> plot(x,y,'rx',x,true) >> x0 = x; >> y0 = y; >> true0 = true; >> h = h/2; >> simpleeuler >> whos Name Size Elements Bytes Density Complex h 1 by Full No i 1 by Full No starty 1 by Full No true 1 by Full No true0 1 by Full No x 1 by Full No x0 1 by Full No y 1 by Full No y0 1 by Full No Grand total is 153 elements using 1224 bytes >> plot(x0,abs(true0-y0),'gx',x,abs(true-y),'yo'); 18

19 function [x] = gausselim(a,b) function [x] = gausselim(a,b) File gausselim.m This subroutine will perform Gaussian elmination on the matrix that you pass to it. i.e., given A and b it can be used to find x, Ax = b To run this file you will need to specify several things: A - matrix for the left hand side. b - vector for the right hand side The routine will return the vector x. ex: [x] = gausselim(a,b) this will perform Gaussian elminiation to find x. N = max(size(a)); Perform Gaussian Elimination for j=2:n, for i=j:n, m = A(i,j-1)/A(j-1,j-1); A(i,:) = A(i,:) - A(j-1,:)*m; b(i) = b(i) - m*b(j-1); 19

20 Perform back substitution x = zeros(n,1); x(n) = b(n)/a(n,n); for j=n-1:-1:1, x(j) = (b(j)-a(j,j+1:n)*x(j+1:n))/a(j,j); >> A = [ ; ; ; ] A = >> b = [ ]' b = >> [x] = gausselim(a,b) x = >> function [x,y] = eulerapprox(startx,h,x,starty,func) file: eulerapprox.m This matlab subroutine will find the approximation to 20

21 a D.E. given by y' = func(x,y) y(startx) = starty To run this file you will first need to specify the following: startx : the starting value for x h : the step size x : the ing value for x starty : the initial value func : routine name to calculate the right hand side of the D.E.. This must be specified as a string. ex: [x,y] = eulerapprox(0,1,1/16,1,'f'); Will return the approximation of a D.E. where x is from 0 to 1 in steps of 1/16. The initial value is 1, and the right hand side is calculated in a subroutine given by f.m. The routine will generate two vectors. The first vector is x which is the grid points starting at x0=0 and have a step size h. The second vector is an approximation to the specified D.E. x = [startx:h:x]; y = 0*x; y(1) = starty; for i=2:max(size(y)), y(i) = y(i-1) + h*feval(func,x(i-1),y(i-1)); function [f] = f(x,y) Evaluation of right hand side of a differential equation. f = 1/y; 21

22 >> help eulerapprox file: eulerapprox.m This matlab subroutine will find the approximation to a D.E. given by y' = func(x,y) y(startx) = starty To run this file you will first need to specify the following: startx : the starting value for x h x : the step size : the ing value for x starty : the initial value func : routine name to calculate the right hand side of the D.E.. This must be specified as a string. ex: [x,y] = eulerapprox(0,1,1/16,1,'f'); Will return the approximation of a D.E. where x is from 0 to 1 in steps of 1/16. The initial value is 1, and the right hand side is calculated in a subroutine given by f.m. The routine will generate two vectors. The first vector is x which is the grid points starting at x0=0 and have a step size h. The second vector is an approximation to the specified D.E. >> [x,y] = eulerapprox(0,1/16,1,1,'f'); >> plot(x,y) 22

23 decision = 3; leftx = 0; rightx = 1; lefty = 1; righty = 1; N= 10; h = (rightx-leftx)/(n-1); x = [leftx:h:rightx]'; A = zeros(n); for i=2:n-1, A(i,i-1:i+1) = [1-2 1]; A = A/h^2; A(1,1) = 1; A(N,N) = 1; b = sin(x); b(1) = lefty; b(n) = righty; if(decision<3) Find and plot the eigen values [e,v] = eig(a); e = diag(e); plot(real(e),imag(e),'rx'); title('eigen Values of the matrix'); elseif(decision>3) Find and plot the eigen values of inv(a) 23

24 [e,v] = eig(inv(a)); e = diag(e); plot(real(e),imag(e),'rx'); title('eigen Values of the inverse of the matrix'); else Solve the system y = A\b; linear = (lefty-righty+sin(leftx)-sin(rightx))/(leftx-rightx); constant = lefty + sin(leftx) - linear*leftx; true = -sin(x) + linear*x + constant; subplot(1,2,1); plot(x,y,'go',x,true,'y'); title('true Solution and Approximation'); xlabel('x'); ylabel('y'); subplot(1,2,2); plot(x,abs(y-true),'cx'); title('error'); xlabel('x'); ylabel(' Error '); 24

ENGG1410-F Tutorial 6

ENGG1410-F Tutorial 6 Jianwen Zhao Department of Computer Science and Engineering The Chinese University of Hong Kong 1/16 Problem 1. Matrix Diagonalization Diagonalize the following matrix: A = [ ] 1 2 4 3 2/16 Solution The

穨matlab教學範例ccc.doc

穨matlab教學範例ccc.doc