Report

Matlab iiI Solving non-linear algebra problems Justin Dawber September 22, 2011 Expectations/PreRequisites • Introduction to MatLab I & II (or equal experience) • • • • • • • MatLab as a calculator Anonymous Functions Function and Script m-files MatLab Basics Element-wise Operations 2-D and 3-D Plotting Code Cells and Publishing Anonymous Function Review • Used in calling functions indirectly • >> Sin = @sin; % The variable ‘Sin’ points to the function ‘sin’ • >> Sin(pi) % Evaluates the sine of pi • (Not the most useful example, more later) • Can be used to create ‘anonymous functions’ • >> myfun = @(x) 1./(x.ˆ3 + 3*x - 5) • >> myfun(3) M-file review • Scripts: No input and no output arguments. Contain a series of commands that may call other scripts and functions • Functions: Accept input and returns output arguments. Usually called program routines and have a special definition syntax. • Code Cells: Defined breaks in code that will help breakdown solution process What defines “non-linear” • Any Equation where the output is not directly proportional to the input (or does not obey the superposition principle) • Simplest Examples = 2 , = , = sin() • Polynomials 3 + 3 2 + 5 = 0 • Single Variable Equations 1.44 = 2 /(1 − )2 • Non-linear Multivariate Equations 2 − − + 2 • The terms make this non-linear • Multivariate linear equations are also possible HUMPS (Built-in) • HUMPS is a built-in Matlab function that has a strong maxima at 0.3 • For those that want to know: humps( x) 1 1 6 2 ( x 0.3) 0.01 ( x 0.9) 2 0.04 Introducing fsolve • Solver for systems of non-linear equations • Requires Optimization Toolbox • Single Dimension System • Use fsolve with an anonymous function • Steps to a solution • Define the Anonymous Function • Plot the function to visualize initial guess • Call fsolve with function and guess • solution = fsolve(@f, initial_guess) Lets See It! • We will do the following: • Review a script m-file in Matlab for HUMPS • • • • • Start with a complete srcipt Review Code cells to break-up the code Plot the function to visualize the guess Iterate through the cells Review the initial guess and solutions • Switch to MatLab An Example: Finding Equilibrium Concentrations • For the reaction H2O + CO <=> CO2 + H2 • Given the Equilibrium constant, K = 1.44 • For an equimolar feed of H2O and CO, compute the equilibrium conversion • Solution • • • • • C_H2O = (C_H20)0 * (1-Xe) C_CO = (C_H20)0 * (1-Xe) C_H2 = (C_H20)0 * Xe C_CO2 = (C_H20)0 * Xe K = C_CO2*C_H2/(C_H20*C_CO) Lets Try It: • We will do the following: • Write a script m-file in Matlab for the Equilibrium Conversion • Start with a skeleton script • Use Code cells to break-up the code • Plot the function to visualize the guess • Review a common syntax problem for element-wise operations • Iterate through the cells • Review the initial guess and solutions • Switch to MatLab 2-Dimensional System of non-linear equations • What do we have in this case? • 2 surfaces • What are we trying to find? • The point in x and y where both surfaces are zero • What is different about this case? • Hard to visualize • Two initial guesses required • Requires a Function-Function m-file • Also know as sub-functions or function in a function The multi-dimensional function m-file • Use sub-functions (function-function) • Primary function – call fsolve • Secondary or sub-function – define the multi-variate system function main clear all; close all; clc; %% We can make some plots to help identify initial guesses x = 0:0.1:2; y=x; [X,Y] = meshgrid(x,y); hold on surf(X,Y,2*X - Y - exp(X) + 2) % first function surf(X,Y,-X + 2*Y - exp(Y) + 2) % second function zlabel('z') view(69,8) %% initial_guesses = [1.3 0.9]; [X,fval,exitflag] = fsolve(@myfunc,initial_guesses) function z = myfunc(X) % set of equations to be solved x = X(1); y = X(2); z = [2*x - y - exp(x) + 2; -x + 2*y - exp(y) + 2]; Lets see it: • We will do the following: • Review a script m-file in Matlab for two stretched HUMPS surfaces • • • • • Start with the complete script Review each of the Code cells Plot the function to visualize the guess Review the initial guess Call fsolve and review solutions • Switch to MatLab An example: 2-Dimensional system • Find a solution to the following system: • 2 − = + 2 • − + 2 = − 2 Lets see it: • We will do the following: • Review a script m-file in Matlab for the preceding example • • • • • Start with the complete script Review each Code cell Plot the function Review the initial guess Call fsolve and review solutions • Switch to MatLab n-Dimensional Systems Example • fsolve also works for the n-dimensional case • No way to graph this case • Must have knowledge of the problem to specify and initial guess • Solve the set of equations: • • • • • • 0 = K1-(yCO*yH2^3)/(yCH4*yH2O)*p^2; 0 = K2 - ((yCO^2)*(yH2^2))/(yCO2*yCH4); 0 = (2*yCH4 + yH2 + yH2O)*n - nH2f; 0 = (0.5*yCO + yCO2 + 0.5*yH2O)*n - nO2f; 0 = (yCH4 + yCO + yCO2)*n - nCf; 0 = yCO + yCO2 + yCH4 + yH2 + yH2O - 1; • Given • • • • K1 = 25.82; K2 = 19.41; p = 1; nCf = 1; nH2f=0.5; nO2f = 0.5; Lets see it: • We will do the following: • Review a script m-file in Matlab for the n-Dim example • • • • • Start with the complete script Review each Code cell Cannot plot the function Postulate on the initial guess Call fsolve and review solutions • Switch to MatLab Thing to consider when using fsolve • No Solutions • Multiple Solutions • Depends on the initial guess • Infinite Solutions – coincidence • The nature of Numerical Solvers – Know your tolerances Lets take a look at one of each: No solution example • Translated HUMPS • Let’s slide the HUMPS graph up 50 • It no longer crosses the X-axis • We can attempt to solve it in the same way • Lets see how fsolve handles it? Lets See it: • We will do the following: • Run the earlier script for the 1-D humps example with the graph translated +50 • • • • Start with the complete script Run the script Confirm the translation of +50 Review the output from fsolve • Switch to MatLab Multiple Solution example • Back to the earlier HUMPS example • Two different guesses yield two different solutions • As you can see, two Zeros. A guess around -.4 will return the lower zero, while a guess near 1.2 will yield the high one. Infinite Solutions Example • Back to a 2-D fsolve example • Solve the system of equations: • • sin(x) + sin(y) - 1.5 -sin(x) - sin(y) + 1.5 Lets See it: • We will review the graph of the two surfaces in the preceding example • View graph from different angles • Call fsolve with multiple initial guesses • Switch to Matlab A little bit about Numerical Solvers Tolerances • Numerical Solvers search for a solution starting from and initial guess • Several search algorithms can be used • Termination criteria • Solver terminates when the value of the function is within a certain range close to 0 • The solver is unaware of the scale of the problem • If you are looking for a value in ml, but the problem is in m3 the solver may stop at ~0.003 m3 … but this is 3 L! • Lets look at an example of this and how to correct it Tolerance Concern Example • Solve this Equation V− (0 − ) 2 =0 • Given • • • • CA0= 2000 mol/m3 V = .010 m3 v = .0005 m3/s k = .00023 m3/mol/s • After plotting we will see solution is near 500 • Guess = 500 • Call fsolve • Guess is exactly right?!? - Something must be wrong • Scaling is the problem Tolerance Concern Example (cont.) • How to Proceed? • Option One – Change the tolerance • Optimset • options = optimset('TolFun',1e-9); • Be careful not to set tolerance too tight (1e-15 = too tight) • Then call fsolve with options • fsolve(f,cguess,options) • Results are now much more accurate Tolerance Concern Example (cont.) • How to proceed? • Option 2 - Scale the input units and guess: • • • • • CA0= 2 mol/L V = 10 L v = .5 L/s k = .23 L/mol/s Guess = .5 • Call fsolve with default tolerances • Results now more accurate Lets see it: • We will do the following: • Review a script m-file in Matlab for the preceding tolerance concern example • • • • Start with the complete script Review each Code cell Iterate through code cells Review solutions using different methods • Switch to MatLab Polynomials in MATLAB • Defining Polynomials • How to get f(x) = Ax2+Bx+C into MatLab • Simple! --- >>f = [A B C] • Finding Roots of a polynomial f • Also Simple --- >>roots(f) • An example: The Van Der Waal equation: • V3 – ((pnb+nRT)/p)V2 + (n2a/p)V – n3ab/p • Coefficients [1 -(pnb+nRT/p) (n2a/p) – n3ab/p] Lets try it: • We will compare the solutions to the above polynomial using fsolve and roots • • • • • • Start with a complete script m-file Define Polynomial as Anonymous Function Define Polynomial as coefficient vector [a b c d] Find solution using roots() and fsolve() roots is an analytical solution • All solutions • fsolve is a numerical solution • Only the Closest Solution • Switch to Matlab Summary • fsolve(@f,initial_guess) – non-linear solver (From the optimization toolbox) • Remember to consider the no/multiple/infinite solution cases • Remember to set your tolerances or scale you problem to ensure accuracy • • • • Provides only the solution closest to the initial guess Quality of initial guess directly related to quality of solution • Intuition about problem is beneficial • Use graphical output to aid in choosing guess Optimset can set various options for solve • • Especially important when using the units package (more on that later) >>help optimset for more info roots() – Solves polynomial defined by their coefficients. • Provides all solutions