Lecture 6

Lecture 6 - Lecture 6. Welcome to the last lecture...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Lecture 6. Welcome to the last lecture on Matlab. This lecture will extend the capabilities of Matlab to solving polynomials as well other calculus applications. This lecture will cover the following topics:  ­ representing and operating on polynomials  ­ numerical derivatives  ­ numerical integrations  ­ function handles  ­ solving ODEs  ­ cell arrays  ­ Representing Polynomials Given the following polynomial p( x ) = x 2 + x ! 6 , we can represent it in Matlab with a coefficient matrix in the form of p = [ 1, 1,  ­6]; where each element in in p corresponds to the coefficient of p(x). For example, the matrix p = [1, 2, 3]; would represent the polynomial p( x ) = x 2 + 2 x + 3 If the polynomial has a degree of n, then the corresponding coefficient matrix will have length a length of n + 1. The last element in the vector is the lowest order term (which is the x 0 term). If you want to find the order of the polynomial, take the length  ­1. All zero coefficients must be marked as a zero in the matrix. This means that if a term does not exist, such as in the equation p( x ) = x 3 ! 2 x ! 5 where the x 2 term is missing, we must represent it in the form of p = [ 1, 0,  ­2,  ­5]; Operating on Polynomials Matlab provides a set of functions to deal with polynomials. The first function we are going to look at is the polyval function. The polyval functions evaluates the result of a polynomial at a given value. The polyval function comes in the form of >>polyval(matrix, value) where matrix is the polynomial and value is what you want x to be. Example: Given the following polynomial, which represents p( x ) = x 2 + 2 x + 3 >> p = [ 1, 2, 3] p = 1 2 3 To evaluate p(1), which is equal to p(1) = 12 + 2 * (1) + 3 , just use >> polyval(p, 1) ans = 6 (Recall the cubic.m script and the polynomial.m scripts that you have written for Assignment 2 and 3. polyval does the exact same thing as those scripts.) We can use the roots function to find the roots of a polynomial. The roots takes the form of >>roots(matrix) where matrix is the coefficient matrix of the polynomial. Example: Given the polynomial x 2 + x ! 6 = 0 >> q = [1, 1,  ­6] q = 1 1  ­6 To find the roots of the polynomial, enter >> roots(q) ans =  ­3 2 Given a set of data, we can find a polynomial to approximately fit the data by using polyfit. The polyfit function determines the coefficients (c0 , c1, c2 , ...cn ) of a polynomial of order n that best fits the data according to y = c0 x n + c1 x n!1 + c2 x n!2 + c3 x n!3 + .. + cn polyfit takes the form of >>polyfit(x,y,n) where x is the matrix for the x values, y is the matrix for the y values, and n is the degree of the polynomial to be returned. Example: Given the following set of data x and y, try to find a polynomial of degree 2 to approximate the equation: >> x = [ 0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]; >> y = [ 2, 0.75, 0,  ­0.25, 0, 0.75, 2]; >> p = polyfit(x,y,2) p = 1.0000  ­3.0000 2.0000 which is y = x 2 ! 3x + 2 We can visualize all of this by plotting the x and y values as well as the polynomial y = x 2 ! 3x + 2 . >> x2 = 0:0.1:3; >> y2 = x2.^2  ­3*x2 + 2; >> plot(x,y,'b ­x',x2,y2,'r ­ ­+') If you have tow polynomials of the same order, you can add them and subtract them using the + and – operators. If they differ in order, you can pad the lower order one with zeros on the left and this will still work. Example: Given the two polynomial, y1 = 2 x 2 + 4 x ! 1 and y2 = 4 x 2 ! 2 x + 2 , we can add and subtract them with the + and the – operators. >> y1 = [2,4, ­1]; >> y2 = [4, ­2,2]; >> y3 = y1 + y2 y3 = 6 2 1 >> y4 = y1  ­ y2 y4 =  ­2 6  ­3 The conv function is equivalent to polynomial multiplication and deconv is equivalent to polynomial division. (Note, deconv returns 2 outputs, which are the answer and the remainder. Type ‘help deconv’ to learn more. Example: Given z1 = x ! 2 and z2 = x + 3 , >> z1 = [1, ­2]; >> z2 = [1,3]; >> z3 = conv(z1,z2) z3 = 1 1  ­6 which is equivalent to z1 * z2 = ( x ! 2 )( x + 3) = x 2 + 3x ! 2 x ! 6 = x 2 + x ! 6 Numerical Derivatives dy Recall that the derivative of a function can be written in the form of , where dy dx and dx are infinitesimal changes in y and in x respectively. To approximate the !y derivative, we can calculate instead. !x Matlab provides a function called diff that will take the difference between adjacent elements. You will get an array with one element shorter the original array. We can calculate diff(y)./diff(x) to approximate the derivative for the data. Example : Let’s try to plot y = x 2 and it’s derivative y' = 2 x : >> x=0:0.1:4; >> y = x.^2; >> x2=0:0.1:3.9; >> y2=diff(y)./diff(x); >> plot(x,y,'g ­',x2,y2,'r ­ ­') >> legend('x^2','derivative of x^2') Note: The reason why we set x2 to be 1 less value than x is because y2 is the result of diff on y and x, which produces 1 less value than x. Similarly, let try to find and plot the derivative of sin(x): >> x = 0:pi/20:2*pi; >> y = sin(x); >> x2 = 0:pi/20:2*pi ­pi/20; >> y2 = diff(y)./diff(x); >> plot(x,y,'g ­',x2,y2,'r ­ ­') >> legend('sin(x)','derivative of sin(x)') Numerical Integration Recall the trapezoid method to find the area under a curve: To find approximate the area under this curve, we can divide the curve into segments, and try to place a trapezoid under each segments and add up the sum. Matlab provides a function called trapz that approximate the integral of a function by approximate using the trapezoid method. The trapz function has the form of >>trapz(x,y) Example: ! Let’s try to approximate the exact value of ! sin( x ) dx . 0 ! We know that " sin( x ) dx = ! cos( x ) ! = ! cos(! ) ! (! cos(0 )) = !(!1) ! (!1) = 2 0 0 >> x = 0:pi/100:pi; >> y = sin(x); >> trapz(x,y) ans = 1.9998 To compute the numerical integral of a function, we can also use the quadrature function, quad. (Quadrature is synonymous with numerical integration. The trapezoid method above is a specific type of quadrature method. The way the quadrature works is by recursively calling a function and evaluates it at certain points and finding the sum. The exact mechanism of how quadrature works is beyond the scope of this lecture.) The quad method takes the form of >>quad(function, a, b) where function is the function handle of the function that you are trying to integrate and a and b indicates the range of the x values. Before we can continue with the discussion of the quad function, we need to learn a little more about function handles. Function Handles We can create a function handle by prepending an @ sign to an EXISTING function. The function can be either built ­in or something you have created and saved in an m ­file. The function handle can be thought of as a reference to the actual function. The most important use of function handles is that we can treat the handle as a variable and pass it into other functions so they can use them too. Example: Given the squaref function from lecture 4: function y = squaref(x) % This function take a value and return the square of it. y= x .^ 2; We can create a function handle for squaref by entering >> s = @squaref; What the function handle does is that provides a means of calling a function indirectly. Now you may use the function handle as though it is the squaref function itself. >> s(10) ans = 100 To use the quad function, we need to pass the function handle to it. Example: Let’s try to integrate the y = x 2 function from 0 to 1 The first thing we need is WRITE the y = x 2 function OURSELVES! Our squaref function does exactly this. So lets just create a function handle for it: >> s = @squaref; Next, let us call the quad function with a = 0 and b = 1 >> quad(s, 0, 1) ans = 0.3333 Let’s calculate the integral by hand and verify the results: 1 111 1 1 ! x 2 dx = 3 x 3 0 = 3 (1)3 " 3 (0) = 3 0 As we can see, the quad function work as intended. Note: Instead of actually creating the function handle separately and then pass it into quad, we could done >>quad(@squaref, 0, 1) and that would have produced the same result. Solving Ordinary Differential Equations (ODEs) dy dy = f ( y, t ) where dt dt is the change in y with respect to time (the derivative of y), and f ( y, t ) is any function of y and time. This is a differential equation because the derivative y is dependent on y itself. dy Here’s an extremely simple differential equation: = ! y . We can solve this ODE by dt separations of variables and then integrate. The resulting equation is y = Ce! t where C is a constant. In order to solve for a specific solution, we must provide an initial condition, such as y(0)=1, so that we can solve for C. In this case, if y(0)=1, C is equal to 1, and our answer is y = e! t . We can solve the ODE in Matlab quite easily using the solver in Matlab. Matlab provide a large set of ODE solvers for us to use, each with it’s own degree of accuracy. The ODE solver we are going to us is ode45. The ODE solver takes the form of >>[t,y] = ode45(odefunction, tspan, y0) dy where odefunction is the handle for the function that gives us the , tspan is the dt range of t values, and y0 is the initial condition. It returns two outputs, t, for the t values, and y, for the y values. dy So let’s try to solve = ! y in Matlab. dt The first thing we need to do is create our ODE fuction, simpleode.m A simple first order differential question has the general form function dy_dt = simpleode(t,y) %Simple ode function dy_dt=-y; Please note, that the order of the inputs (t,y) matters. If the function is instead declared as simpleode(y,t), you will get an incorrect result. The next step is to call the ode45 function. Remember that the ode45 produces 2 outputs. >>[t,y] = ode45(@simpleode, [0,5],[1]); Finally, let’s plot our y vs t graph: >> plot(t,y) Note: This graph should make sense, as the solution to the ODE is y = e! t , which is an exponentially decaying function. Recall from Professor Vallancourt’s lecture on solving ODEs: Given a mass on a spring, we can express the force, F, on the spring with Newton’s d2x Law F = m * 2 and Hooke’s Law F =  ­k * x. Setting the two equations together, we dt 2 dx d2x d2x k get m * 2 = !k * x . Solve for 2 and we get 2 = ! * x . However, the ODE dt dt dt m dx solver requires our function to give the output in the form of = f (t, x ) , so we need dt d 2 x d dx = ( ) . We can rewrite our dt 2 dt dt dy problem into two equations by introducing a new intermediate variable, , where dt dy dy dx k = x . So now we have two equations, = x and = ! * y . Notice that if I take dt dt dt m d2x k the derivative of both sides with respect to dt, we will get back 2 = ! * x . dt m $! x $ 0 d ! y $ !1 So now with the two equations, we can write the matrix # & = # &# & 'k / m%" y % dt " x % "0 to rewrite the equation a little. Recall that function statevars = spring_oscillator(t, vars) % m = 2; k = 5000; y = vars(1); x = vars(2); % dy_dt =x; dx_dt = -k*y/m; % statevars = [dy_dt; dx_dt]; Now let’s call the ode45 function and plot the result: >> [t,x]=ode45(@spring_oscillator, [0,1],[0.1;0]); >> plot(t,x(:,1)) Note: We are setting the initial condition to be that y(0)=0.1 and x(0) = 0 (You can read more on solving and plotting ODEs here http://laser.ceb.cam.ac.uk/wiki/images/e/e5/NumMeth_Handout_7.pdf) Cell Arrays In addition to matrices, Matlab also supports something called a cell array, which is a collection of variables. The biggest difference between a cell array and a matrix, is that a cell array can store matrices of different sizes, where as a matrix can only store rows/columns of the same size. We can also mix matrix with strings and store them in the cell array. To create a cell array, use the curly brackets { } instead of the regular bracket. >>cell = { item1, item2, item2…} To access the items in the cell array, we use the curly brackets to index them. >>cell{1} Example: >> x = {1:10, 'hello', true} x = [1x10 double] 'hello' [1] >> x{1} ans = 1 2 3 4 5 6 7 8 9 10 >> x{2} ans = hello >> x{3} ans = 1 ...
View Full Document

This note was uploaded on 11/03/2011 for the course MATH 1090 taught by Professor Greenwood during the Spring '08 term at MIT.

Ask a homework question - tutors are online