non-linear-algebra

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

similar documents