Matlab - PowerPoint Presentation

Slide1

Matlab iiI

Solving non-linear algebra problems

Justin Dawber

September 22, 2011Slide2

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 PublishingSlide3

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

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 processSlide5

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

 Slide6

HUMPS (Built-in)

HUMPS is a built-in Matlab function that has a strong maxima at 0.3

For those that want to know:Slide7

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 guessCall fsolve with function and guesssolution = fsolve(@f, initial_guess)Slide8

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 guessIterate through the cellsReview the initial guess and solutionsSwitch to MatLabSlide9

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 conversionSolutionC_H2O = (C_H20)0 * (1-Xe)C_CO = (C_H20)0 * (1-Xe)C_H2 = (C_H20)0 * XeC_CO2 = (C_H20)0 * XeK = C_CO2*C_H2/(C_H20*C_CO)Slide10

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 codePlot the function to visualize the guessReview a common syntax problem for element-wise operationsIterate through the cellsReview the initial guess and solutionsSwitch to MatLabSlide11

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 functionSlide12

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

mainclear all; close all; clc;%% We can make some plots to help identify initial guessesx = 0:0.1:2; y=x;[X,Y] = meshgrid(x,y);hold onsurf(X,Y,2*X - Y - exp(X) + 2) % first functionsurf(X,Y,-X + 2*Y - exp(Y) + 2) % second functionzlabel('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];Slide13

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 guessReview the initial guessCall fsolve and review solutionsSwitch to MatLabSlide14

An example: 2-Dimensional system

Find a solution to the following system:

 Slide15

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 functionReview the initial guessCall fsolve and review solutionsSwitch to MatLabSlide16

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;GivenK1 = 25.82;K2 = 19.41;p = 1;nCf = 1; nH2f=0.5; nO2f = 0.5;Slide17

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 functionPostulate on the initial guessCall fsolve and review solutionsSwitch to MatLabSlide18

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:Slide19

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?Slide20

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 MatLabSlide21

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

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.5Slide23

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

MatlabSlide24

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 m

3

the solver may stop at ~0.003 m

3 … but this is 3 L!Lets look at an example of this and how to correct itSlide25

Tolerance Concern Example

Solve this Equation

= 0

Given

C

A0

= 2000

mol

/m

3V = .010 m3v = .0005 m3/sk = .00023 m3/mol/sAfter plotting we will see solution is near 500Guess = 500Call fsolveGuess is exactly right?!? - Something must be wrongScaling is the problem Slide26

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 optionsfsolve(f,cguess,options)Results are now much more accurateSlide27

Tolerance Concern Example (cont.)

How to proceed?

Option 2 - Scale the input units and guess:

C

A0

=

2

mol

/L

V = 10 Lv = .5 L/sk = .23 L/mol/sGuess = .5Call fsolve with default tolerancesResults now more accurateSlide28

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 methodsSwitch to MatLabSlide29

Polynomials in MATLAB

Defining Polynomials

How to get f(x) = Ax

2

+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/pCoefficients [1 -(pnb+nRT/p) (n2a/p) – n3ab/p]Slide30

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 solutionAll solutionsfsolve is a numerical solutionOnly the Closest SolutionSwitch to MatlabSlide31

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 solutionIntuition about problem is beneficialUse graphical output to aid in choosing guessOptimset can set various options for solve>>help optimset for more inforoots() – Solves polynomial defined by their coefficients. Provides all solutions