MIT6_094IAP10_lec03

6.094 Introduction to programming in MATLAB Danilo Šćepanović IAP 2008 Lecture 3 : Solving Equations and Curve Fitting Homework 2 Recap • How long did it take? • Using min with matrices: » a=[3 7 5;1 9 10; 30 -1 2]; » b=min(a); % returns the min of each column » m=min(b); % returns min of entire a matrix » m=min(min(a)); % same as above » m=min(a(:)); % makes a a vector, then gets min • Common mistake: » [m,n]=find(min(a)); % think about what happens • How to make and run a function: save the file, then call it from the command window like any other function. No need to 'compile' or make it official in any other way Outline (1) Linear Algebra (2) Polynomials (3) Optimization (4) Differentiation/Integration (5) Differential Equations Systems of Linear Equations • Given a system of linear equations ¾ x+2y-3z=5 ¾ -3x-y+z=-8 ¾ x-y+z=0 • Construct matrices so the system is described by Ax=b » A=[1 2 -3;-3 -1 1;1 -1 1]; » b=[5;-8;0]; • And solve with a single line of code! » x=A\b; ¾ x is a 3x1 vector containing the values of x, y, and z • The \ will work with square or rectangular systems. • Gives least squares solution for rectangular systems. Solution depends on whether the system is over or underdetermined. MATLAB makes linear algebra fun! More Linear Algebra • Given a matrix » mat=[1 2 -3;-3 -1 1;1 -1 1]; • Calculate the rank of a matrix » r=rank(mat); ¾ the number of linearly independent rows or columns • Calculate the determinant » d=det(mat); ¾mat must be square ¾ if determinant is nonzero, matrix is invertible • Get the matrix inverse » E=inv(mat); ¾ if an equation is of the form A*x=b with A a square matrix, x=A\b is the same as x=inv(A)*b Matrix Decompositions • MATLAB has built-in matrix decomposition methods • The most common ones are » [V,D]=eig(X) ¾ Eigenvalue decomposition » [U,S,V]=svd(X) ¾ Singular value decomposition » [Q,R]=qr(X) ¾QR decomposition Exercise: Linear Algebra • Solve the following systems of equations: ¾ System 1: ¾ System 2: 4 34 3 2 x y x y + = − + = 2 2 4 3 3 4 2 x y x y x y − = − + = + = Exercise: Linear Algebra • Solve the following systems of equations: ¾ System 1: ¾ System 2: 4 34 3 2 x y x y + = − + = 2 2 4 3 3 4 2 x y x y x y − = − + = + = » A=[1 4;-3 1]; » b=[34;2]; » rank(A) » x=inv(A)*b; » A=[2 -2;-1 1;3 4]; » b=[4;3;2]; » rank(A) ¾ rectangular matrix » x1=A\b; ¾ gives least squares solution » error=abs(A*x1-b) Outline (1) Linear Algebra (2) Polynomials (3) Optimization (4) Differentiation/Integration (5) Differential Equations Polynomials • Many functions can be well described by a high-order polynomial • MATLAB represents a polynomials by a vector of coefficients ¾ if vector P describes a polynomial ax3+bx2+cx+d • P=[1 0 -2] represents the polynomial x2-2 • P=[2 0 0 0] represents the polynomial 2x3 P(1) P(2) P(3) P(4) Polynomial Operations • P is a vector of length N+1 describing an N-th order polynomial • To get the roots of a polynomial » r=roots(P) ¾ r is a vector of length N • Can also get the polynomial from the roots » P=poly(r) ¾ r is a vector length N • To evaluate a polynomial at a point » y0=polyval(P,x0) ¾ x0 is a single value; y0 is a single value • To evaluate a polynomial at many points » y=polyval(P,x) ¾ x is a vector; y is a vector of the same size Polynomial Fitting • MATLAB makes it very easy to fit polynomials to data • Given data vectors X=[-1 0 2] and Y=[0 -1 3] » p2=polyfit(X,Y,2); ¾ finds the best second order polynomial that fits the points (-1,0),(0,-1), and (2,3) ¾ see help polyfit for more information » plot(X,Y,’o’, ‘MarkerSize’, 10); » hold on; » x = -3:.01:3; » plot(x,polyval(p2,x), ‘r--’); Exercise: Polynomial Fitting • Evaluate for x=-4:0.1:4. • Add random noise to these samples. Use randn. Plot the noisy signal with . markers • Fit a 2nd degree polynomial to the noisy data • Plot the fitted polynomial on the same plot, using the same x values and a red line 2y x= Exercise: Polynomial Fitting • Evaluate for x=-4:0.1:4. » x=-4:0.1:4; » y=x.^2; • Add random noise to these samples. Use randn. Plot the noisy signal with . markers » y=y+randn(size(y)); » plot(x,y,’.’); • Fit a 2nd degree polynomial to the noisy data » p=polyfit(x,y,2); • Plot the fitted polynomial on the same plot, using the same x values and a red line » hold on; » plot(x,polyval(p,x),’r’) 2y x= Outline (1) Linear Algebra (2) Polynomials (3) Optimization (4) Differentiation/Integration (5) Differential Equations Nonlinear Root Finding • Many real-world problems require us to solve f(x)=0 • Can use fzero to calculate roots for any arbitrary function • fzero needs a function passed to it. • We will see this more and more as we delve into solving equations. • Make a separate function file » x=fzero('myfun',1) » x=fzero(@myfun,1) ¾ 1 specifies a point close to where you think the root is Courtesy of The MathWorks, Inc. Used with permission. Minimizing a Function • fminbnd: minimizing a function over a bounded interval » x=fminbnd('myfun',-1,2); ¾myfun takes a scalar input and returns a scalar output ¾myfun(x) will be the minimum of myfun for -1≤x ≤ 2 • fminsearch: unconstrained interval » x=fminsearch('myfun',.5) ¾ finds the local minimum of myfun starting at x=0.5 Anonymous Functions • You do not have to make a separate function file » x=fzero(@myfun,1) ¾What if myfun is really simple? • Instead, you can make an anonymous function » x=fzero(@(x)(cos(exp(x))+x^2-1), 1 ); » x=fminbnd(@(x) (cos(exp(x))+x^2-1),-1,2); input function to evaluate Optimization Toolbox • If you are familiar with optimization methods, use the optimization toolbox • Useful for larger, more structured optimization problems • Sample functions (see help for more info) » linprog ¾ linear programming using interior point methods » quadprog ¾ quadratic programming solver » fmincon ¾ constrained nonlinear optimization Exercise: Min-Finding • Find the minimum of the function over the range –π to π. Use fminbnd. • Plot the function on this range to check that this is the minimum. ( ) ( ) ( )cos 4 sin 10 xf x x x e−= Exercise: Min-Finding • Find the minimum of the function over the range –π to π. Use fminbnd. • Plot the function on this range to check that this is the minimum. • Make the following function: » function y=myFun(x) » y=cos(4*x).*sin(10*x).*exp(-abs(x)); • Find the minimum in the command window: » x0=fminbnd('myFun',-pi,pi); • Plot to check if it's right » figure; x=-pi:.01:pi; plot(x,myFun(x)); ( ) ( ) ( )cos 4 sin 10 xf x x x e−= Outline (1) Linear Algebra (2) Polynomials (3) Optimization (4) Differentiation/Integration (5) Differential Equations Numerical Differentiation • MATLAB can 'differentiate' numerically » x=0:0.01:2*pi; » y=sin(x); » dydx=diff(y)./diff(x); ¾ diff computes the first difference • Can also operate on matrices » mat=[1 3 5;4 8 6]; » dm=diff(mat,1,2) ¾ first difference of mat along the 2nd dimension, dm=[2 2;4 -2] ¾ see help for more details ¾ The opposite of diff is the cumulative sum cumsum • 2D gradient » [dx,dy]=gradient(mat); 0 100 200 300 400 500 600 700 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Numerical Integration • MATLAB contains common integration methods • Adaptive Simpson's quadrature (input is a function) » q=quad('myFun',0,10); ¾ q is the integral of the function myFun from 0 to 10 » q2=quad(@(x) sin(x)*x,0,pi) ¾ q2 is the integral of sin(x)*x from 0 to pi • Trapezoidal rule (input is a vector) » x=0:0.01:pi; » z=trapz(x,sin(x)); ¾ z is the integral of sin(x) from 0 to pi » z2=trapz(x,sqrt(exp(x))./x) ¾ z2 is the integral of from 0 to pixe x Outline (1) Linear Algebra (2) Polynomials (3) Optimization (4) Differentiation/Integration (5) Differential Equations ODE Solvers: Method • Given a differential equation, the solution can be found by integration: ¾ Evaluate the derivative at a point and approximate by straight line ¾ Errors accumulate! ¾ Variable timestep can decrease the number of iterations ODE Solvers: MATLAB • MATLAB contains implementations of common ODE solvers • Using the correct ODE solver can save you lots of time and give more accurate results » ode23 ¾ Low-order solver. Use when integrating over small intervals or when accuracy is less important than speed » ode45 ¾High order (Runge-Kutta) solver. High accuracy and reasonable speed. Most commonly used. » ode15s ¾ Stiff ODE solver (Gear's algorithm), use when the diff eq's have time constants that vary by orders of magnitude ODE Solvers: Standard Syntax • To use standard options and variable time step » [t,y]=ode45('myODE',[0,10],[1;0]) • Inputs: ¾ODE function name (or anonymous function). This function takes inputs (t,y), and returns dy/dt ¾ Time interval: 2-element vector specifying initial and final time ¾ Initial conditions: column vector with an initial condition for each ODE. This is the first input to the ODE function • Outputs: ¾ t contains the time points ¾ y contains the corresponding values of the integrated variables. ODE integrator: 23, 45, 15s ODE function Time range Initial conditions ODE Function • The ODE function must return the value of the derivative at a given time and function value • Example: chemical reaction ¾ Two equations ¾ODE file: – y has [A;B] – dydt has [dA/dt;dB/dt] A B 10 50 10 50 10 50 dA A B dt dB A B dt = − + = − Courtesy of The MathWorks, Inc. Used with permission. ODE Function: viewing results • To solve and plot the ODEs on the previous slide: » [t,y]=ode45('chem',[0 0.5],[0 1]); ¾ assumes that only chemical B exists initially » plot(t,y(:,1),'k','LineWidth',1.5); » hold on; » plot(t,y(:,2),'r','LineWidth',1.5); » legend('A','B'); » xlabel('Time (s)'); » ylabel('Amount of chemical (g)'); » title('Chem reaction'); ODE Function: viewing results • The code on the previous slide produces this figure 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (s) A m o u n t o f c h e m i c a l ( g ) Chem reaction A B Higher Order Equations • Must make into a system of first-order equations to use ODE solvers • Nonlinear is OK! • Pendulum example: ( ) ( ) ( ) 0g sin L g sin L let g sin L x dx dt θ θ θ θ θ γ γ θ θ γ θ γ + = = − = = − ⎡ ⎤ = ⎢ ⎥⎣ ⎦ ⎡ ⎤ = ⎢ ⎥⎣ ⎦ && && & & v v & & Courtesy of The MathWorks, Inc. Used with permission. Plotting the Output • We can solve for the position and velocity of the pendulum: » [t,x]=ode45('pendulum',[0 10],[0.9*pi 0]); ¾ assume pendulum is almost horizontal » plot(t,x(:,1)); » hold on; » plot(t,x(:,2),'r'); » legend('Position','Velocity'); 0 1 2 3 4 5 6 7 8 9 10 -8 -6 -4 -2 0 2 4 6 8 Position Velocity Position in terms of angle (rad) Velocity (m/s) Plotting the Output • Or we can plot in the phase plane: » plot(x(:,1),x(:,2)); » xlabel('Position'); » yLabel('Velocity'); • The phase plane is just a plot of one variable versus the other: -3 -2 -1 0 1 2 3 -8 -6 -4 -2 0 2 4 6 8 Position V e l o c i t y Velocity is greatest when theta=0 Velocity=0 when theta is the greatest ODE Solvers: Custom Options • MATLAB's ODE solvers use a variable timestep • Sometimes a fixed timestep is desirable » [t,y]=ode45('chem',[0:0.001:0.5],[0 1]); ¾ Specify the timestep by giving a vector of times ¾ The function value will be returned at the specified points ¾ Fixed timestep is usually slower because function values are interpolated to give values at the desired timepoints • You can customize the error tolerances using odeset » options=odeset('RelTol',1e-6,'AbsTol',1e-10); » [t,y]=ode45('chem',[0 0.5],[0 1],options); ¾ This guarantees that the error at each step is less than RelTol times the value at that step, and less than AbsTol ¾Decreasing error tolerance can considerably slow the solver ¾ See doc odeset for a list of options you can customize Exercise: ODE • Use ode45 to solve for on the range t=[0 10], with initial condition and • Plot the result. 10dy dt t y= − ( )y t ( )0 10y = Exercise: ODE • Use ode45 to solve for on the range t=[0 10], with initial condition and • Plot the result. • Make the following function » function dydt=odefun(t,y) » dydt=-t*y/10; • Integrate the ODE function and plot the result » [t,y]=ode45(‘odefun’,[0 10],10); • Alternatively, use an anonymous function » [t,y]=ode45(@(t,y) –t*y/10,[0 10],10); • Plot the result » plot(t,y);xlabel('Time');ylabel('y(t)'); 10dy dt t y= − ( )y t ( )0 10y = Exercise: ODE • The integrated function looks like this: 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 Time y ( t ) Function y(t), integrated by ode45 End of Lecture 3 (1) Linear Algebra (2) Polynomials (3) Optimization (4) Differentiation/Integration (5) Differential Equations We're almost done! MIT OpenCourseWare http://ocw.mit.edu 6.094 Introduction to MATLAB® January (IAP) 2010 For information about citing these materials or our Terms of Use, visit: http://ocw.mit.edu/terms. 6.094�Introduction to programming in MATLAB Homework 2 Recap Outline Systems of Linear Equations More Linear Algebra Matrix Decompositions Exercise: Linear Algebra Exercise: Linear Algebra Outline Polynomials Polynomial Operations Polynomial Fitting Exercise: Polynomial Fitting Exercise: Polynomial Fitting Outline Nonlinear Root Finding Minimizing a Function Anonymous Functions Optimization Toolbox Exercise: Min-Finding Exercise: Min-Finding Outline Numerical Differentiation Numerical Integration Outline ODE Solvers: Method ODE Solvers: MATLAB ODE Solvers: Standard Syntax ODE Function ODE Function: viewing results ODE Function: viewing results Higher Order Equations Plotting the Output Plotting the Output ODE Solvers: Custom Options Exercise: ODE Exercise: ODE Exercise: ODE End of Lecture 3 << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles true /AutoRotatePages /All /Binding /Left /CalGrayProfile (None) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.4 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJDFFile false /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.1000 /ColorConversionStrategy /sRGB /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams false /MaxSubsetPct 100 /Optimize true /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo false /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Remove /UsePrologue false /ColorSettingsFile () /AlwaysEmbed [ true ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 150 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 150 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.50000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /ColorImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 150 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /GrayImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 300 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects true /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /False /Description << /JPN /DEU /FRA