Using mesh function (which plots a surface defined by a series of points in 3d)

ID: 3766369 • Letter: U

Question

Using mesh function (which plots a surface defined by a series of points in 3d) combined with the meshgrid function (which creates 2d grid defined by two vectors), show the surface defined by the following equation—known as a saddle function: f(x, y) = (x – 3)2 – (y – 2)2 for 2 < x < 4 and 1 < y < 3 (increment: 0.2).

Using mesh function (which plots a surface defined by a series of points in 3d) combined with the meshgrid function (which creates 2d grid defined by two vectors), show the surface defined by the following equation—known as a saddle function: f(x, y) = (x – 3)2 – (y – 2)2 for 2 < x < 4 and 1 < y < 3 (increment: 0.2).

Explanation / Answer

To assign a value to a variable in MATLAB simply type the name of the variable, followed by the assignment operator, =, followed by the value.

>> x=7

x =

Note that variable names in MATLAB are case sensitive, so X and x are not equal. We can perform all of the usual operations with x.

>> x^2-3*x+2

ans =

>> log(x)

ans =

1.94591014905531

>> sin(x)

ans =

0.65698659871879

New variables may be defined using existing variables.

>> y=8*x

y =

This, however, does not imply any permanent relationship between x and y. If we change x, the value of y does not change.

>> x=x+5

x =

>> y

y =

The command who returns a list of all variables in the current workspace, while whos returns the same list with more detailed information about each variable.

>> who

Your variables are:

ans x y

>> whos

Name Size Bytes Class

ans 1x1 8 double array

x 1x1 8 double array

y 1x1 8 double array

Grand total is 3 elements using 24 bytes

Notice that the size of each variable is 1×1. All variables in MATLAB are matrices. Scalars such as x and y are just 1×1 matrices. We will explain how to enter matrices in the next section. To clear one or more variables from the workspace, type clear followed by the names of the variables. Typing just clear clears all variables.

>> clear

>> who

>> x

??? Undefined function or variable 'x'.

4 Matrices and Vectors

To enter a matrix in MATLAB, use square brackets and separate entries within a row by spaces and separate rows using semicolons.

>> A=[2 1 -1 8; 1 0 8 -3; 7 1 2 4]

A =

2 1 -1 8

1 0 8 -3

7 1 2 4

Often we do not want MATLAB to display a response, especially when dealing with very large matrices. To suppress the output, place a semicolon at the end of the line. Typing

>> B=[2 0 -3; -1 1 3];

will still define the variable B containing a 2×3 matrix, but MATLAB will not echo anything.

>> whos

Name Size Bytes Class

A 3x4 96 double array

B 2x3 48 double array

v 3x1 24 double array

Grand total is 21 elements using 168 bytes

To view the contents of the variable B, just type its name.

>> B

B =

2 0 -3

-1 1 3

Vectors (column vectors) are simply matrices with a single column.

>> v = [ 2; 3; -4]

v =

-4

A row vector is a matrix with a single row.

>> w=[3 -2 5 11]

w =

3 -2 5 11

It is often necessary to define vectors with evenly spaced entries. In MATLAB, the colon (:) provides a shorthand for creating such vectors.

>> 2:5

ans =

2 3 4 5

Typing j:i:k defines a row vector with increment i starting at j and ending at k.

>> 3:2:9

ans =

3 5 7 9

Recall that the transpose of a matrix A is the matrix AT whose entry in row i column j is the same as the entry in row j column i of A. In MATLAB, A' represents the transpose of the matrix A.

>> A=[5 -2 9; 11 7 8]

A =

5 -2 9

11 7 8

>> A'

ans =

5 11

-2 7

9 8

To define regularly spaced column vectors, we can take the transpose of a regularly spaced row vector.

>> [1:3:10]'

ans =

The entry in row i, column j of a matrix A is A(i,j).

>> A=[3 -2 7 8; 4 3 2 1; 10 15 -2 9]

A =

3 -2 7 8

4 3 2 1

10 15 -2 9

>> A(3,2)

ans =

It is also possible to view multiple entries within any row or column. For instance, the second and fourth entries in the third row are accessed as follows.

>> A(3,[2 4])

ans =

15 9

Row i of A is A(i,:) and column j of A is A(:,j).

>> A(3,:)

ans =

10 15 -2 9

>> A(:,3)

ans =

-2

Next we display the first, second and fourth columns.

>> A(:,[1 2 4])

ans =

3 -2 8

4 3 1

10 15 9

The entries of a vector (row or column) may be accessed using a single index.

>> w=[7; 13; 11]

w =

>> w(2)

ans =

Matrices with the same number of rows may be concatenated horizontally, and matrices with the same number of columns may be concatenated vertically.

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

A =

1 2 3

4 5 6

>> B=[7 8; 9 10]

B =

7 8

9 10

>> [A B]

ans =

1 2 3 7 8

4 5 6 9 10

>> C=[7 8 9]

C =

7 8 9

>> [A;C]

ans =

1 2 3

4 5 6

7 8 9

To remove rows or columns from a matrix, simply redefine them to be empty matrices.

>> A=[ 4 7 2 1 3; 8 7 12 -2 5; 11 1 14 -2 0]

A =

4 7 2 1 3

8 7 12 -2 5

11 1 14 -2 0

>> A(2,:)=[]

A =

4 7 2 1 3

11 1 14 -2 0

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

A =

7 1 3

1 -2 0

5 Dot Products and Cross Products

The function dot computes the dot product of two vectors in Rn .

>> v=[7; 23; 15; 2], w=[5; -2; 1; -8]

v =

w =

-2

-8

>> dot(v,w)

ans =

-12

Note that the dot product is symmetric.

>> dot(w,v)

ans =

-12

Recall that the length of a vector v is ||v||={v·v}.

>> vlength=sqrt(dot(v,v))

vlength =

28.4077

The length of a vector can also be found directly using the norm function.

>> norm(v)

ans =

28.4077

Also recall that if q is the angle between two vectors v and w then v·w=||v||||w||cosq. Solving for the angle we have q = arccos((v·w)/||v||||w||).

>> theta=acos(dot(v,w)/(norm(v)*norm(w)))

theta =

1.6144

>> theta*180/pi

ans =

92.4971

So the angle between v and w is about 92.5.

The function cross computes the cross product of two vectors in R3.

>> v=[3; 2; 1], w=[4; 15; 1]

v =

w =

>> x=cross(v,w)

x =

-13

Cross products are anti-symmetric. That is, w×v=-v×w.

>> cross(w,v)

ans =

-1

-37

The cross product v×w is orthogonal to both v and w. We can verify this by taking its dot product with both v and w. Recall that two vectors are orthogonal if and only if their dot product equals zero.

>> dot(x,v)

ans =

>> dot(x,w)

ans =

6 Basic Matrix Operations

Addition (and subtraction) of matrices of the same dimensions is performed componentwise.

>> A=[5 -1 2; 3 4 7]

A =

5 -1 2

3 4 7

>> B=[2 2 1; 5 0 3]

B =

2 2 1

5 0 3

>> A+B

ans =

7 1 3

8 4 10

Note that only matrices with the same dimensions may be summed.

>> C=[3 1; 6 4]

C =

3 1

6 4

>> A+C

??? Error using ==> +

Matrix dimensions must agree.

Scalar multiplication (and division by nonzero scalars) is also performed componentwise.

>> 2*A

ans =

10 -2 4

6 8 14

The * in scalar multiplication is not optional.

>> 2A

??? 2

Missing operator, comma, or semi-colon.

Vector addition and scalar multiplication are performed in the same way.

>> v=[3; 5], w=[-2; 7]

v =

w =

-2

>> 10*v-5*w

ans =

The matrix product A*B is defined when A is m×n and B is n×k. That is, the number of columns of A must equal the number of rows of A. In this case the product A*B is an m×k matrix.

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

A =

3 1 7 2

6 -3 4 2

9 4 -1 -2

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

B =

1 2

3 4

5 6

7 8

>> A*B

ans =

55 68

31 40

2 12

The entry in row i, column j of A*B is the dot product of row i of A with column j of B.

>> dot(A(2,:),B(:,1))

ans =

MATLAB produces an error message if the inner matrix dimensions do not agree.

>> B*A

??? Error using ==> *

Inner matrix dimensions must agree.

The exception to this rule is when one of the matrices is a 1×1 matrix, i.e. a scalar. In this case the product is interpreted as scalar multiplication.

>> C=[2]

C =

>> A*C

ans =

6 2 14 4

12 -6 8 4

18 8 -2 -4

Matrix-vector products are special cases of matrix products.

>> A=[13 -11 21; 16 9 10], v=[19; -7; 15]

A =

13 -11 21

16 9 10

v =

-7

>> A*v

ans =

639

391

As above, the entries of A*v are the dot products of the rows of A with v.

>> [dot(A(1,:),v); dot(A(2,:),v)]

ans =

639

391

Also, the product A*v equals the linear combination of the columns of A whose coefficients are the components of v.

>> v(1)*A(:,1)+v(2)*A(:,2)+v(3)*A(:,3)

ans =

639

391

A square matrix may be multiplied by itself, and in this case it makes sense to take powers of the matrix. For instance, A^6 equals A*A*A*A*A*A.

>> A=[0 1; 1 1]

A =

0 1

1 1

>> A*A

ans =

1 1

1 2

>> A^6

ans =

5 8

8 13

7 Reduced Row Echelon Form

The rref command is used to compute the reduced row echelon form of a matrix.

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

A =

1 2 3 4 5 6

1 2 4 8 16 32

2 4 2 4 2 4

1 2 1 2 1 2

>> rref(A)

ans =

1 2 0 0 -4 -8

0 0 1 0 -1 -6

0 0 0 1 3 8

0 0 0 0 0 0

Any problem that can be phrased in terms of a system of linear equations can thus be solved using rref. For instance, consider the following set of vectors.

>> v1=[1; 1; 1]; v2=[1; -2; 1]; v3=[1; 2; 3]; v4=[2; 3; 4]; v5=[1; -1; -3];

Suppose we want to write v5 as a linear combination of the other vectors. First we define the augmented matrix for the resulting system of equations x1v1+x2v2+x3v3+x4v4=v5.

>> A=[v1 v2 v3 v4 v5]

A =

1 1 1 2 1

1 -2 2 3 -1

1 1 3 4 -3

The reduced row echelon form then yields the solution(s).

>> rref(A)

ans =

1 0 0 1 3

0 1 0 0 0

0 0 1 1 -2

The variable x4 is a free variable, so let's choose x4=0. Then x1=3, x2=0 and x3=-2, and therefore v5=3v1-2v3. Let's check.

3*v1-2*v3

ans =

-1

-3

Any other choice of x4 will result in a different way of writing v5 as a linear combination of the other vectors. Now suppose we want to write v2 as a linear combination of v1, v3, v4 and v5.

>> A=[v1 v3 v4 v5 v2]

A =

1 1 2 1 1

1 2 3 -1 -2

1 3 4 -3 1

>> rref(A)

ans =

1 0 1 3 0

0 1 1 -2 0

0 0 0 0 1

This time the third equation reduces to 0=1, so there are no solutions, and therefore v2 is not a linear combination of the other vectors.

Another way to use rref is in the form [R,p]=rref(A), which defines R to be the reduced row echelon form of A and p to be the vector listing the columns of R which contain pivots.

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

A =

1 2 1 5

1 2 2 6

1 2 3 7

1 2 4 8

>> [R,p]=rref(A)

R =

1 2 0 4

0 0 1 1

0 0 0 0

p =

1 3

Extracting the pivot columns of A gives a basis for the column space of A.

>> A(:,p)

ans =

1 1

1 2

1 3

1 4

With the reduced row echelon form of A in hand we could easily find a basis for the null space of A. The command null does this for us.

>> null(A,'r')

ans =

-2 -4

1 0

0 -1

0 1

The 'r' stands for rational and tells MATLAB to find the null space from the reduced row echelon form of A. Without the 'r', MATLAB finds an orthonormal basis for the null space - that is, a basis consisting of mutually orthogonal unit vectors.

>> N=null(A)

N =

-0.9608 0

0.1601 -0.8165

-0.1601 -0.4082

0.1601 0.4082

To verify that this is an orthonormal basis we form the product N'*N. Since the rows of N' are the columns of N, the entry in row i, column j of N'*N is the dot product of column i and column j of N.

>> N'*N

ans =

1.0000 0

0 1.0000

The ones on the diagonal indicate that the columns of N are unit vectors, and the zeros indicate that they are orthogonal to each other.

We can also use rref to find the inverse of an invertible matrix. Recall that a matrix A is invertible if and only if rref(A)=In, the n×n identity matrix, and that in this case rref[A | In]=[In | A-1]. For example, consider the following matrix.

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

A =

1 1 1

1 2 3

1 3 6

In MATLAB, the n×n identity matrix In is given by eye(n). Let's augment A by eye(3) and compute its reduced row echelon form.

>> B=[A eye(3)]

B =

1 1 1 1 0 0

1 2 3 0 1 0

1 3 6 0 0 1

>> rref(B)

ans =

1 0 0 3 -3 1

0 1 0 -3 5 -2

0 0 1 1 -2 1

>> Ainv=[3 -3 1; -3 5 -2; 1 -2 1]

Ainv =

3 -3 1

-3 5 -2

1 -2 1

is the inverse of A. To verify this, recall that A-1A=AA-1=In.

>> Ainv*A

ans =

1 0 0

0 1 0

0 0 1

>> A*Ainv

ans =

1 0 0

0 1 0

0 0 1

8 Rank

To compute the rank of a matrix we use the rank command. We briefly recall some of the important facts regarding the rank of a matrix.

For example,

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

A =

1 2 1 4

2 3 1 3

3 2 1 2

4 3 1 1

>> rank(A)

ans =

Thus the column space of A has dimension 3 and the columns of A above are linearly dependent. Using rank we can determine which columns of A form a basis for its column space.

>> rank(A(:,[1 2 3]))

ans =

>> rank(A(:,[1 2 4]))

ans =

>> rank(A(:,[1 3 4]))

ans =

>> rank(A(:,[2 3 4]))

ans =

Thus any choice of three columns forms a basis for the column space of A except columns one, three and four. Next we use rank as a test for invertibility.

>> A=[11 -21 3; 8 2 1; 16 -12 5]

A =

11 -21 3

8 2 1

16 -12 5

>> rank(A)

ans =

>> B=[3 4 5;6 7 8;9 10 11]

B =

3 4 5

6 7 8

9 10 11

>> rank(B)

ans =

The matrix A is invertible, but B is not.

Using rank it is possible to determine whether or not a given vector is in the column space of a matrix.

>> A=[5 8 -4; 3 19 11; -6 6 0; 12 4 1]

A =

5 8 -4

3 19 11

-6 6 0

12 4 1

>> v1=[21; 16; -7; 33],v2=[30; 7; 30; -16]

v1 =

-7

v2 =

-16

>> rank(A)

ans =

So the column space of A is three-dimensional. Now consider the augmented matrix [A v]. If v is in the column space of A, then [A v] has rank 3. Otherwise [A v] has rank 4.

>> rank([A v1])

ans =

>> rank([A v2])

ans =

Thus v2 is in the column space of A, but v1 is not.

9 Inverses

The inverse of an invertible matrix A can be found by using either A^(-1) or inv(A).

>> A=[2 1 1; 1 2 2; 2 1 2]

A =

2 1 1

1 2 2

2 1 2

>> Ainv=inv(A)

Ainv =

2/3 -1/3 0

2/3 2/3 -1

-1 0 1

Let's verify the result.

>> A*Ainv

ans =

1 0 0

0 1 0

0 0 1

>> Ainv*A

ans =

1 0 0

0 1 0

0 0 1

MATLAB gives a warning message if the matrix is singular (not invertible).

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

B =

1 2 3

4 5 6

7 8 9

>> inv(B)

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 2.055969e-018.

ans =

1.0e+016 *

-0.4504 0.9007 -0.4504

0.9007 -1.8014 0.9007

-0.4504 0.9007 -0.4504

Here MATLAB actually returns an inverse. To see that B is really singular, compute its rank.

>> rank(B)

ans =

Since the rank of B is less than 3, B is singular. Also recall that a matrix is invertible if and only if its determinant is nonzero.

>> det(B)

ans =

Unfortunately, there are invertible matrices which MATLAB regards as singular.

>> format long

>> C=[1.00000000000001 1; 1 .99999999999999]

C =

1.00000000000001 1.00000000000000

1.00000000000000 0.99999999999999

>> inv(C)

Warning: Matrix is singular to working precision.

ans =

Inf Inf

>> rank(C)

ans =

>> det(C)

ans =

>> rref(C)

ans =

1.00000000000000 0.99999999999999

0 0

According to MATLAB, the matrix C has rank 1 and determinant zero, and is therefore singular. However, if we let e = 0.00000000000001, then det(C)=(1+e)(1-e)-1=-e2 0, so C is invertible. The problem is that in format long MATLAB is accurate to 15 places and therefore recognizes 1+e and 1-e as being different from 1, but recognizes

(1+e)(1-e)=1-e2 = 0.9999999999999999999999999999

as 1. This example should serve as a warning to not blindly accept everything that MATLAB tells us.

To solve an equation of the form Ax=b where A is invertible, we simply multiply by the inverse of A to get x=A-1b.

>> A=[11 7 -6 8; 3 -1 12 15; 1 1 1 7; -4 6 1 8]

A =

11 7 -6 8

3 -1 12 15

1 1 1 7

-4 6 1 8

>> b=[10; -23; -13; 4]

b =

-23

-13

>> format rat

>> x=inv(A)*b

x =

-3

Let's verify this result.

>> A*x

ans =

-23

-13

10 Eigenvectors and Eigenvalues

The eig command is used to find the eigenvalues of a square matrix.

>> A=[ 3 1 1; 1 3 1; 1 1 3]

A =

3 1 1

1 3 1

1 1 3

>> eig(A)

ans =

2.0000

5.0000

If A is diagonalizable, the eig command can also be used to find an eigenbasis, along with the diagonal matrix to which it is similar.

>> [Q,D]=eig(A)

Q =

-0.8164 -0.0137 0.5774

0.3963 0.7139 0.5774

0.4201 -0.7001 0.5774

D =

2.0000 0 0

0 2.0000 0

0 0 5.0000

Here, the columns of the matrix Q form an eigenbasis of A, and Q-1AQ=D. Let's check.

>> inv(Q)*A*Q

ans =

2.0000 0 0.0000

0.0000 2.0000 0.0000

-0.0000 0.0000 5.0000

The matrix Q actually contains an orthonormal basis of eigenvectors.

>> Q'*Q

ans =

1.0000 0.0000 -0.0000

0.0000 1.0000 -0.0000

-0.0000 -0.0000 1.0000

If we just wanted to find the eigenvectors in the ''usual'' way, we would use the null command.

>> C1=null(A-2*eye(3),'r')

C1 =

-1 -1

1 0

0 1

>> C2=null(A-5*eye(3),'r')

C2 =

Let's now check that these three vectors form an eigenbasis.

>> C=[C1 C2]

C =

-1 -1 1

1 0 1

0 1 1

>> inv(C)*A*C

ans =

2.0000 -0.0000 0.0000

0 2.0000 -0.0000

0 0 5.0000

Great! It worked.

17 Solving Algebraic Equations

The solve command is used to find solutions of equations involving symbolic expressions.

>> solve('sin(x)+x=5')

ans =

5.6175550052726989176213921571114

In expressions with more than one variable, we can solve for one or more of the variables in terms of the others. Here we find the roots of the quadratic ax2+bx+c in x in terms of a, b and c. By default solvesets the given expression equal to zero if an equation is not given.

>> solve('a*x^2+b*x+c','x')

ans =

[ 1/2/a*(-b+(b^2-4*a*c)^(1/2))]

[ 1/2/a*(-b-(b^2-4*a*c)^(1/2))]

Systems of equations can also be handled by solve.

>> S=solve('x+y+z=1','x+2*y-z=3')

S =

x: [1x1 sym]

y: [1x1 sym]

The variable S contains the solution, which consists of x and y in terms of z.

>> S.x

ans =

-3*z-1

>> S.y

ans =

2*z+2

Now let's find the points of intersection of the circles x2+y2=4 and (x-1)2+(y-1)2=1.

>> S=solve('x^2+y^2=4','(x-1)^2+(y-1)^2=1')

S =

x: [2x1 sym]

y: [2x1 sym]

>> [S.x S.y]

ans =

[ 5/4-1/4*7^(1/2), 5/4+1/4*7^(1/2)]

[ 5/4+1/4*7^(1/2), 5/4-1/4*7^(1/2)]

The points of intersection are therefore ((5-7)/4,(5+7)/4) and ((5+7)/4,(5-7)/4).

18 Derivatives

Differentiation of a symbolic expression is performed by means of the function diff. For instance, let's find the derivative of f(x)=sin(ex).

>> syms x

>> f=sin(exp(x))

f =

sin(exp(x))

>> diff(f)

ans =

cos(exp(x))*exp(x)

The nth derivative of f is diff(f,n).

>> diff(f,2)

ans =

-sin(exp(x))*exp(x)^2+cos(exp(x))*exp(x)

To compute the partial derivative of an expression with respect to some variable we specify that variable as an additional argument in diff. Let f(x,y)=x3y4+ysinx.

>> syms x y

>> f=x^3*y^4+y*sin(x)

f =

x^3*y^4+y*sin(x)

First we compute f/x.

>> diff(f,x)

ans =

3*x^2*y^4+y*cos(x)

Next we compute f/y.

>> diff(f,y)

ans =

4*x^3*y^3+sin(x)

Finally we compute 3f/x3.

>> diff(f,x,3)

ans =

6*y^4-y*cos(x)

The Jacobian matrix of a function f:RnRm can be found directly using the jacobian function. For example, let f:R2R3 be defined by f(x,y)=(sin(xy),x2+y2,3x-2y).

>> f=[sin(x*y); x^2+y^2; 3*x-2*y]

f =

[ sin(y*x)]

[ x^2+y^2]

[ 3*x-2*y]

>> Jf=jacobian(f)

Jf =

[ cos(y*x)*y, cos(y*x)*x]

[ 2*x, 2*y]

[ 3, -2]

In the case of a linear transformation, the Jacobian is quite simple.

>> A=[11 -3 14 7;5 7 9 2;8 12 -6 3]

A =

11 -3 14 7

5 7 9 2

8 12 -6 3

>> syms x1 x2 x3 x4

>> x=[x1;x2;x3;x4]

x =

[ x1]

[ x2]

[ x3]

[ x4]

>> T=A*x

T =

[ 11*x1-3*x2+14*x3+7*x4]

[ 5*x1+7*x2+9*x3+2*x4]

[ 8*x1+12*x2-6*x3+3*x4]

Now let's find the Jacobian of T.

>> JT=jacobian(T)

JT =

[ 11, -3, 14, 7]

[ 5, 7, 9, 2]

[ 8, 12, -6, 3]

The Jacobian of T is precisely A.

Next suppose f:RnR is a scalar valued function. Then its Jacobian is just its gradient. (Well, almost. Strictly speaking, they are the transpose of one another since the Jacobian is a row vector and the gradient is a column vector.) For example, let f(x,y)=(4x2-1)e-x2-y2.

>> syms x y real

>> f=(4*x^2-1)*exp(-x^2-y^2)

f =

(4*x^2-1)*exp(-x^2-y^2)

>> gradf=jacobian(f)

gradf =

[ 8*x*exp(-x^2-y^2)-2*(4*x^2-1)*x*exp(-x^2-y^2), -2*(4*x^2-1)*y*exp(-x^2-y^2)]

Next we use solve to find the critical points of f.

>> S=solve(gradf(1),gradf(2));

>> [S.x S.y]

ans =

[ 0, 0]

[ 1/2*5^(1/2), 0]

[ -1/2*5^(1/2), 0]

Thus the critical points are (0,0), (5/2,0) and (-5/2,0).

The Hessian of a scalar valued function f:RnR is the n×n matrix of second order partial derivatives of f. In MATLAB we can obtain the Hessian of f by computing the Jacobian of the Jacobian of f. Consider once again the function f(x,y)=(4x2-1)e-x2-y2.

>> syms x y real

>> Hf=jacobian(jacobian(f));

>> Hf=simple(Hf)

Hf =

[2*exp(-x^2-y^2)*(2*x+1)*(2*x-1)*(2*x^2-5), 4*x*y*exp(-x^2-y^2)*(-5+4*x^2)]

[4*x*y*exp(-x^2-y^2)*(-5+4*x^2), 2*exp(-x^2-y^2)*(-1+2*y^2)*(2*x+1)*(2*x-1)]

We can now use the Second Derivative Test to determine the type of each critical point of f found above.

>> subs(Hf,{x,y},{0,0})

ans =

10 0

0 2

>> subs(Hf,{x,y},{1/2*5^(1/2),0})

ans =

-5.7301 0

0 -2.2920

>> subs(Hf,{x,y},{-1/2*5^(1/2),0})

ans =

-5.7301 0

0 -2.2920

Thus f has a local minimum at (0,0) and local maxima at the other two critical points. Evaluating f at the critical points gives the maximum and minimum values of f.

>> subs(f,{x,y},{0,0})

ans =

-1

>> subs(f,{x,y},{'1/2*5^(1/2)',0})

ans =

4*exp(-5/4)

>> subs(f,{x,y},{'-1/2*5^(1/2)',0})

ans =

4*exp(-5/4)

Thus the minimum value of f is f(0,0)=-1 and the maximum value is f(5/2,0)=f(-5/2,0)=4e-5/4. The graph of f is shown below.

As our final example, we solve a Lagrange multiplier problem. For f(x,y)=xy(1+y) let's find the maximum and minimum of f on the unit circle x2+y2=1. First we enter the function f and the constraint functiong(x,y)=x2+y2-1.

>> syms x y mu

>> f=x*y*(1+y)

f =

x*y*(1+y)

>> g=x^2+y^2-1

g =

x^2+y^2-1

Next we solve the Lagrange multiplier equations f(x,y)-mg(x,y)=0 and constraint equation g(x,y)=0 for x, y and m.

>> L=jacobian(f)-mu*jacobian(g)

L = [ y*(1+y)-2*mu*x, x*(1+y)+x*y-2*mu*y]

>> S=solve(L(1),L(2),g)

S =

mu: [5x1 sym]

x: [5x1 sym]

y: [5x1 sym]

Next let's view the critical points found. We can ignore m now.

>> [S.x S.y]

ans =

[ 1/6*(22-2*13^(1/2))^(1/2), 1/6+1/6*13^(1/2)]

[ -1/6*(22-2*13^(1/2))^(1/2), 1/6+1/6*13^(1/2)]

[ 1/6*(22+2*13^(1/2))^(1/2), 1/6-1/6*13^(1/2)]

[ -1/6*(22+2*13^(1/2))^(1/2), 1/6-1/6*13^(1/2)]

[ 0, -1]

Next we need to evaluate f at each of these points.

>> values=simple(subs(f,{x,y},{S.x,S.y}))

values =

[ 1/216*(22-2*13^(1/2))^(1/2)*(1+13^(1/2))*(7+13^(1/2))]

[ -1/216*(22-2*13^(1/2))^(1/2)*(1+13^(1/2))*(7+13^(1/2))]

[ 1/216*(22+2*13^(1/2))^(1/2)*(-1+13^(1/2))*(-7+13^(1/2))]

[ -1/216*(22+2*13^(1/2))^(1/2)*(-1+13^(1/2))*(-7+13^(1/2))]

[ 0]

Finally we convert these into decimal expressions to identify the maximum and minimum. This is done using the double command.

>> double(values)

ans =

0.8696

-0.8696

-0.2213

0.2213

Thus the maximum of f is about 0.8696 and the minimum is about -0.8696.

19 M-Files

% <-- The % is for comments.

% MATLAB ignores everything following it on the same line.

% test.m

A=[1 2; 3 4]

B=[5 6; 7 8]

A+B

After creating a file containing the text above, we save it with the filename test.m. In MATLAB we access this file by typing test. Note that, in order for this to work we must be sure that the current directory in MATLAB contains the file. The Unix commands ls (list contents of current directory), cd (change directory) and pwd (list name of current directory) all work within MATLAB.

>> test

A =

1 2

3 4

B =

5 6

7 8

ans =

6 8

10 12

MATLAB executes the commands in the M-file in order, as if we had typed them within MATLAB. We can also do loops within a script.

% sumsquares.m

% sums the first n squares up to n=10

s=0;

for n=1:10

s=s+n^2

end

Here is what this looks like in MATLAB.

>> sumsquares

s =

140

s =

204

s =

285

s =

385

Here's one last example to try. It shows an animation of the curve r(t)=(2tcost/(t+1),2tsint/(t+1)) being plotted. See the Plotting Curves section for details on the commands used.

% curve.m

% Shows animation of a parametric curve being plotted.

hold on

for T=0:.1:4*pi

t=[T T+.1];

plot(2*t.*cos(t)./(t+1),2*t.*sin(t)./(t+1))

axis equal

axis([-2 2 -2 2])

axis off

pause(.01)

end

(1+e)(1-e)=1-e2 = 0.9999999999999999999999999999

Navigate

Using matlab: 3.2 Simpson\'s Rule you can access the interpolation polynomial gi

Using methods on the first picture, can someone help me write a java program to

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.

Using mesh function (which plots a surface defined by a series of points in 3d)

Question

Explanation / Answer

Related Questions

Navigate