Matlab iiI

1. MatlabiiI Solving non-linear algebra problems Justin Dawber September 22, 2011

2. 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

3. 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)

4. 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

5. 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 • Polynomials • Single Variable Equations • Non-linear Multivariate Equations • The terms make this non-linear • Multivariate linear equations are also possible

6. HUMPS (Built-in) • HUMPS is a built-in Matlab function that has a strong maxima at 0.3 • For those that want to know:

7. 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)

8. 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

9. 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)

10. 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

11. 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

12. 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); holdon 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) functionz = 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];

13. 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

14. An example: 2-Dimensional system • Find a solution to the following system:

15. 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

16. 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;

17. 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

18. 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:

19. 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?

20. 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

21. 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.

22. 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

23. 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

24. 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

25. Tolerance Concern Example • Solve this Equation = 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

26. 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

27. 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

28. 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

29. 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]

30. 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

31. 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 • Especially important when using the units package (more on that later) • 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 • >>help optimset for more info • roots() – Solves polynomial defined by their coefficients. • Provides all solutions