Transcript Solution of Nonlinear Equation Systems

Solution of Nonlinear Equation Systems
F(x)  0
Daniel Baur
ETH Zurich, Institut für Chemie- und Bioingenieurwissenschaften
ETH Hönggerberg / HCI F128 – Zürich
E-Mail: [email protected]
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / Nonlinear Equations
1
Zero of Nonlinear Equation Systems
 Problem definition:
 Find the solution of F(x) = 0, where both F and x are vector
(or matrix) valued; Look for the solution either
 In an interval xlb < x < xub
 Or in the uncertainty interval [a, b]
 Types of algorithm available:
1. Substitution algorithms
2. Methods based on function approximation
 Assumptions:
 At least one zero exists in the defined interval
 We are looking for one zero, not all of them
Daniel Baur / Numerical Methods for Chemical Engineers / Nonlinear Equations
2
Function linearization
 f1  x , y  
f x  
0
 f2  x, y 
 Let us again consider Taylor expansion
f1 ( x , y )  f1 ( x 0 , y 0 ) 
f 2 ( x , y )  f 2 ( x0 , y0 ) 
 f1
x
( x  x0 ) 
( x0 , y 0 )
f2
x
( x  x0 ) 
( x0 , y 0 )
 f1
y
( y  y0 )  0
( x0 , y 0 )
f2
y
( y  y0 )  0
( x0 , y 0 )
 In matrix form, this reads
  f1 /  x

f2 / x
 f1 /  y 

f 2 / y  x , y
 0 0
 x  x0 
 f1 ( x 0 , y 0 ) 






y

y
f
(
x
,
y
)
0 

 2 0 0 
 Which is equivalent to
J (x0 ) x   f (x0 )

Newton Method!
Daniel Baur / Numerical Methods for Chemical Engineers / Nonlinear Equations
3
The Jacobian matrix
 Consider a continuous, differentiable function
J i ,k 
fi
xk
  f1
 x
 1
 f2

J   x1



 fn
  x1
 f1
x2
f 2
 x1
f n
x2
 f1 
xn 

f 2 

xn



f n 
 x n 
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
f :
n

 f 1  x1 , x 2 , x 3 ,

f x ,x ,x ,
 2 1 2 3


 f n  x1 , x 2 , x 3 ,
n
, xn  

, xn 



, x n  
4
The Multidimensional Newton-Algorithm
1. Define a starting point x
2. Compute –f(x)
3. Compute J(x)

Either analytically or numerically
4. Solve the linear system of equations J(x)*Δx = -f(x)
for Δx, then set x = x + Δx
 Iterate 2 through 4 until some stopping criteria are fulfilled
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
5
Improvements on the Algorithm
 It might not be necessary to compute the Jacobian at every
iteration, the algorithm can also work nicely if the Jacobian
is only computed every few iterations
 When doing this, instead of solving the linear system of
equations, it might be more efficient to compute inv(J)
when recomputing J, or factorize it in a convenient manner
 A numerical Jacobian can sometimes lead to better results
than an analytical one
 Another option is to adapt the step size, so before
accepting a Δx, one checks if norm(f(x+Δx)) is in fact
smaller than norm(f(x)), and if not, a step of size λ*Δx
with λ < 1 is tried instead
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
6
How does Matlab do it? Nonlinear Systems
 fsolve is the general solver for non-linear equation
systems; It can use different algorithms, namely
 Trust-region algorithms: Instead of just going in the direction that the
Newton step suggests, a check is performed if xk+1 is really better
than xk. This is done by approximating the behavior of the function
around xk (the trust-region). An optimization is then performed to find
an optimal step.
 Trust-region-dogleg: The so-called dogleg method is used to solve
the optimization problem, which constructs the step as a convex
combination of the Cauchy step (steepest descent) and a GaussNewton step.
 Trust-region-reflective: The optimization problem is solved via a
conjugate gradient method in a 2-D subspace.
 The Levenberg-Marquardt method, which finds the step by solving
J(x
) J ( xk )  k I  d k   J ( xk ) f ( xk )
k
T
T
Daniel Baur / Numerical Methods for Chemical Engineers / Nonlinear Equations
7
Matlab Syntax Hints
 fsolve uses the syntax
 x = fsolve(nl_fun, x0, options);
 where nl_fun is a function taking as input x, returning as output the
function values f at x, both can be vectors or matrices
 x0 is an initial guess
 options can be set by options = optimset(...) to choose
algorithms, give analytical Jacobians, etc.; see doc fsolve for
details
Daniel Baur / Numerical Methods for Chemical Engineers / Nonlinear Equations
8
Exercise: CSTR with multiple reactions
 Consider a CSTR where two reactions take place
in
c A,B , Q
in
1
A+B  
C
k
2
C+B  
D
k
ci , Q
out
 We assume the following
 V = const., i.e. Qin = Qout = const.
 Constant temperature and density of the reaction mixture
 Constant feed, i.e. cA,in = const., cB,in = const.
Daniel Baur / Numerical Methods for Chemical Engineers / Nonlinear Equations
9
CSTR mass balance
(Accumulation) = Flow (In – Out) + Reaction
A:
d
dt
B:
d
dt
C:
d
dt
D:
d
dt
 V c A   Q  c Ain
 cA

 V c B   Q  c Bin
 c B  V   k1c A c B  k 2 c B c C 
V   k1c A c B 
V cC   Q  0  cC 
V  k1c A c B  k 2 c B c C 
V c D   Q  0  c D 
V  k 2 c B c C 
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
10
CSTR steady state
d
x   cA
 0
dt
cB
cC
cD 
A:
0   x1  x1 
B:
0   x 2  x 2      k 1 x1 x 2  k 2 x 2 x 3 
C:
0   x3
   k 1 x1 x 2  k 2 x 2 x 3 
D:
0   x4
  k 2 x 2 x 3 
in
    k 1 x1 x 2 
in
where   V / Q
is the residence time
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
11
Assignment 1
1. Show for the 1D-Newton method that it converges
quadratically (locally, order p = 2)
 Remember that for smooth, p+1 times differentiable fixed point
formulas if the following holds: If
 0     x      x  

 p

x  0


 p 1
x
then Φ converges locally with order p
 The Newton method in fixed point form reads
x k 1    x k   x k 
 xk 
f   xk 
f
 Calculate analytically Φ’ and Φ’’, then substitue a fixed point for xk
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
12
Assignment 2
1. Write down the analytical Jacobian matrix for the steady
state CSTR shown on slide 11.
2. Implement the basic Newton method as shown on slide 5.
 Use a function of the form
function [x, info] = newtonMethod(f, J, x0, tol);
where f is a function handle to the function you want to solve, J is a
function handle that returns the Jacobian matrix, x0 is an initial
guess and tol is a vector of tolerances.
 As in with the secand method, use a while loop to find the solution.
 Suggest stopping criteria and failure checks. When can the Newton
method fail in general?
 Use left division «\» to solve the linear system at every iteration (do
not use inv(J)!)
 Let info be a struct you can use to return additional information, like
reason of termination and number of steps needed.
Daniel Baur / Numerical Methods for Chemical Engineers / Nonlinear Equations
13
Assignment 2 (continued)
3. Use your Newton algorithm to solve the steady state
CSTR numerically.
 Create two function files; one that calculates the CSTR equations (1)
as functions of x, and one that calculates the analytical Jacobian as
a function of x.
 Use xIn1 = 1; xIn2 = 1.5; k1 = 0.5; k2 = 10; and τ = 5
 What is the total conversion of A to D?
 Compare your result to what fsolve() finds. Try different starting
guesses. Can you find more than one solution?
4. Find online the function jacobianest.
 Modify your Newton algorithm so that it uses jacobianest to
approximate the Jacobian if the input J is empty (use isempty(J)
to check). To provide an empty input, use [] in the call.
 How many steps are required with the analytical Jacobian compared
to the numerical Jacobian? Which algorithm takes longer?
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
14