Jacobian Matrices Catch Up
undefined
J
a
c
o
b
i
a
n
M
a
t
r
i
c
e
s
C
a
t
c
h
U
p
1
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
S
S
a
a
r
r
a
a
h
h
D
D
u
u
c
c
l
l
o
o
s
s
I
I
v
v
e
e
t
t
i
i
c
c
h
h
E
E
T
T
H
H
Z
Z
u
u
r
r
i
i
c
c
h
h
,
,
I
I
n
n
s
s
t
t
i
i
t
t
u
u
t
t
f
f
ü
ü
r
r
C
C
h
h
e
e
m
m
i
i
e
e
-
-
u
u
n
n
d
d
B
B
i
i
o
o
i
i
n
n
g
g
e
e
n
n
i
i
e
e
u
u
r
r
w
w
i
i
s
s
s
s
e
e
n
n
s
s
c
c
h
h
a
a
f
f
t
t
e
e
n
n
E
E
T
T
H
H
H
H
ö
ö
n
n
g
g
g
g
e
e
r
r
b
b
e
e
r
r
g
g
/
/
H
H
C
C
I
I
F
F
1
1
0
0
6
6
–
–
Z
Z
ü
ü
r
r
i
i
c
c
h
h
E
E
-
-
M
M
a
a
i
i
l
l
:
:
s
s
a
a
r
r
a
a
h
h
.
.
d
d
u
u
c
c
l
l
o
o
s
s
@
@
c
c
h
h
e
e
m
m
.
.
e
e
t
t
h
h
z
z
.
.
c
c
h
h
h
h
t
t
t
t
p
p
s
s
:
:
/
/
/
/
s
s
h
h
i
i
h
h
l
l
a
a
b
b
.
.
e
e
t
t
h
h
z
z
.
.
c
c
h
h
/
/
e
e
d
d
u
u
c
c
a
a
t
t
i
i
o
o
n
n
/
/
S
S
n
n
m
m
/
/
n
n
u
u
m
m
e
e
r
r
i
i
c
c
s
s
.
.
h
h
t
t
m
m
l
l
Z
e
r
o
o
f
N
o
n
l
i
n
e
a
r
E
q
u
a
t
i
o
n
S
y
s
t
e
m
s
P
r
o
b
l
e
m
d
e
f
i
n
i
t
i
o
n
:
F
i
n
d
t
h
e
s
o
l
u
t
i
o
n
o
f
F
(
x
)
=
0
,
w
h
e
r
e
b
o
t
h
F
a
n
d
x
a
r
e
v
e
c
t
o
r
(
o
r
m
a
t
r
i
x
)
v
a
l
u
e
d
;
L
o
o
k
f
o
r
t
h
e
s
o
l
u
t
i
o
n
e
i
t
h
e
r
I
n
a
n
i
n
t
e
r
v
a
l
x
l
b
<
x
<
x
u
b
T
y
p
e
s
o
f
a
l
g
o
r
i
t
h
m
a
v
a
i
l
a
b
l
e
:
1.
Substitution algorithms
2.
Methods based on function approximation
A
s
s
u
m
p
t
i
o
n
s
:
A
t
l
e
a
s
t
o
n
e
z
e
r
o
e
x
i
s
t
s
i
n
t
h
e
d
e
f
i
n
e
d
i
n
t
e
r
v
a
l
W
e
a
r
e
l
o
o
k
i
n
g
f
o
r
o
n
e
z
e
r
o
,
n
o
t
a
l
l
o
f
t
h
e
m
2
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
F
u
n
c
t
i
o
n
l
i
n
e
a
r
i
z
a
t
i
o
n
Let us again consider Taylor expansion
In matrix form, this reads
Which is equivalent to
3
Newton method
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
T
h
e
J
a
c
o
b
i
a
n
m
a
t
r
i
x
Consider a continuous, differentiable function
4
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
T
h
e
M
u
l
t
i
d
i
m
e
n
s
i
o
n
a
l
N
e
w
t
o
n
-
A
l
g
o
r
i
t
h
m
1.
Define a starting point
x
x
2.
Compute
–f(x)
–f(x)
3.
Compute
J
J
(x)
(x)
Either analytically or numerically
4.
Solve the linear system of equations
J
J
(x)*
(x)*
Δ
Δ
x = -f(x)
x = -f(x)
for
Δ
Δ
x
x
, then set
x = x +
x = x +
Δ
Δ
x
x
Iterate 2 through 4 until some stopping criteria are fulfilled
5
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
O
r
d
i
n
a
r
y
D
i
f
f
e
r
e
n
t
i
a
l
E
q
u
a
t
i
o
n
s
(
O
D
E
s
)
6
Sarah Duclos Ivetich / Numerical Methods for Chemical Engineers / Explicit ODE Solvers
S
S
a
a
r
r
a
a
h
h
D
D
u
u
c
c
l
l
o
o
s
s
I
I
v
v
e
e
t
t
i
i
c
c
h
h
E
E
T
T
H
H
Z
Z
u
u
r
r
i
i
c
c
h
h
,
,
I
I
n
n
s
s
t
t
i
i
t
t
u
u
t
t
f
f
ü
ü
r
r
C
C
h
h
e
e
m
m
i
i
e
e
-
-
u
u
n
n
d
d
B
B
i
i
o
o
i
i
n
n
g
g
e
e
n
n
i
i
e
e
u
u
r
r
w
w
i
i
s
s
s
s
e
e
n
n
s
s
c
c
h
h
a
a
f
f
t
t
e
e
n
n
E
E
T
T
H
H
H
H
ö
ö
n
n
g
g
g
g
e
e
r
r
b
b
e
e
r
r
g
g
/
/
H
H
C
C
I
I
F
F
1
1
0
0
6
6
–
–
Z
Z
ü
ü
r
r
i
i
c
c
h
h
E
E
-
-
M
M
a
a
i
i
l
l
:
:
s
s
a
a
r
r
a
a
h
h
.
.
d
d
u
u
c
c
l
l
o
o
s
s
@
@
c
c
h
h
e
e
m
m
.
.
e
e
t
t
h
h
z
z
.
.
c
c
h
h
h
h
t
t
t
t
p
p
s
s
:
:
/
/
/
/
s
s
h
h
i
i
h
h
l
l
a
a
b
b
.
.
e
e
t
t
h
h
z
z
.
.
c
c
h
h
/
/
e
e
d
d
u
u
c
c
a
a
t
t
i
i
o
o
n
n
/
/
S
S
n
n
m
m
/
/
n
n
u
u
m
m
e
e
r
r
i
i
c
c
s
s
.
.
h
h
t
t
m
m
l
l
P
r
o
b
l
e
m
D
e
f
i
n
i
t
i
o
n
We are going to look at initial value problems of
explicit
ODE systems, which means that they can be cast in the
following form
7
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
E
x
a
m
p
l
e
:
A
1
-
D
F
i
r
s
t
O
r
d
e
r
O
D
E
Consider the following initial value problem
This ODE and its solution follow a general form
(Johann Bernoulli, Basel, 17
th
– 18
th
century)
8
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
S
o
m
e
N
o
m
e
n
c
l
a
t
u
r
e
On the next slides, the following notation will be used
9
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
N
u
m
e
r
i
c
a
l
S
o
l
u
t
i
o
n
o
f
a
1
-
D
F
i
r
s
t
O
r
d
e
r
O
D
E
Since we know the first derivative, Taylor expansion
suggests itself as a first approach
This is known as the forward / explicit Euler method,
named after Leonhard Euler (Basel, 18
th
century)
10
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
D
e
f
i
n
i
t
i
o
n
o
f
a
n
i
m
p
l
i
c
i
t
a
l
g
o
r
i
t
h
m
In an implicit algorithm, the function value at the next step
appears on the right hand side of the step equation:
A simple example is the backward / implicit Euler method:
A system of equations has to be solved at every iteration
step;
Why even use implicit algorithms?
11
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
H
o
w
d
o
e
s
M
a
t
l
a
b
d
o
i
t
?
N
o
n
-
S
t
i
f
f
S
o
l
v
e
r
s
ode45
ode45
: Most general solver, best to try first; Uses an
explicit Runge-Kutta pair (Dormand-Prince, p = 4 and 5)
ode23
ode23
: More efficient than
ode45
ode45
at crude tolerances,
better with moderate stiffness; Uses an explicit Runge-
Kutta pair (Bogacki-Shampine, p = 2 and 3).
ode113
ode113
: Better at stringent tolerances and when the ODE
function file is very expensive to evaluate; Uses a variable
order Adams-Bashforth-Moulton method
12
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
H
o
w
d
o
e
s
M
a
t
l
a
b
d
o
i
t
?
S
t
i
f
f
S
o
l
v
e
r
s
ode15s
ode15s
: Use this if
ode45
ode45
fails or is very inefficient and
you suspect that the problem is stiff; Uses multi-step
numerical differentiation formulas (approximates the
derivative in the next point by extrapolation from the
previous points)
ode23s
ode23s
may be more efficient than
ode15s
ode15s
at crude
tolerances and for some problems (especially if there are
largely different time-scales / dynamics involved); Uses a
modified Rosenbrock formula of order 2 (one-step method)
13
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
M
a
t
l
a
b
S
y
n
t
a
x
H
i
n
t
s
All ODE-solvers use the syntax
[T, Y] = ode45(ode_fun, tspan, y0, …);
[T, Y] = ode45(ode_fun, tspan, y0, …);
o
d
e
_
f
u
n
i
s
a
f
u
n
c
t
i
o
n
h
a
n
d
l
e
t
a
k
i
n
g
t
w
o
i
n
p
u
t
s
,
a
s
c
a
l
a
r
t
(
c
u
r
r
e
n
t
t
i
m
e
p
o
i
n
t
)
a
n
d
a
v
e
c
t
o
r
y
(
c
u
r
r
e
n
t
f
u
n
c
t
i
o
n
v
a
l
u
e
s
)
,
a
n
d
r
e
t
u
r
n
i
n
g
a
s
o
u
t
p
u
t
a
v
e
c
t
o
r
o
f
t
h
e
d
e
r
i
v
a
t
i
v
e
s
d
y
/
d
t
i
n
t
h
e
s
e
p
o
i
n
t
s
t
s
p
a
n
i
s
e
i
t
h
e
r
a
t
w
o
c
o
m
p
o
n
e
n
t
v
e
c
t
o
r
[
t
s
t
a
r
t
,
t
e
n
d
]
o
r
a
v
e
c
t
o
r
o
f
t
i
m
e
p
o
i
n
t
s
w
h
e
r
e
t
h
e
s
o
l
u
t
i
o
n
i
s
n
e
e
d
e
d
[
t
s
t
a
r
t
,
t
1
,
t
2
,
.
.
.
]
y
0
a
r
e
t
h
e
v
a
l
u
e
s
o
f
y
a
t
t
s
t
a
r
t
U
s
e
p
a
r
a
m
e
t
r
i
z
i
n
g
f
u
n
c
t
i
o
n
s
t
o
p
a
s
s
m
o
r
e
a
r
g
u
m
e
n
t
s
t
o
o
d
e
_
f
u
n
i
f
n
e
e
d
e
d
f
f
_
_
n
n
e
e
w
w
=
=
@
@
(
(
t
t
,
,
y
y
)
)
o
o
d
d
e
e
_
_
f
f
u
u
n
n
(
(
t
t
,
,
y
y
,
,
k
k
,
,
p
p
)
)
;
;
This can be done directly in the solver call
[T, Y] = ode45(@(t,y)ode_fun(t, y, k, p), tspan,
[T, Y] = ode45(@(t,y)ode_fun(t, y, k, p), tspan,
y0);
y0);
14
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
A
s
s
i
g
n
m
e
n
t
1
1.
Write down the analytical Jacobian matrix for the steady
state CSTR;
2.
Implement the basic Newton method as shown.
Use a function of the form
function
function
[x, info] = newtonMethod(f, J, x0, tol);
[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(
inv(
J
J
)
)
!)
Let
info
info
be a struct you can use to return additional information, like
reason of termination and number of steps needed.
15
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
A
s
s
i
g
n
m
e
n
t
1
(
c
o
n
t
i
n
u
e
d
)
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 x
In
1
= 1; x
In
2
= 1.5; k
1
= 0.5; k
2
= 10; and
τ
= 5
What is the total conversion of A to D?
Compare your result to what
fsolve()
fsolve()
finds. Try different starting
guesses. Can you find more than one solution?
4.
Find online the function
jacobianest
jacobianest
.
M
o
d
i
f
y
y
o
u
r
N
e
w
t
o
n
a
l
g
o
r
i
t
h
m
s
o
t
h
a
t
i
t
u
s
e
s
j
j
a
a
c
c
o
o
b
b
i
i
a
a
n
n
e
e
s
s
t
t
t
o
a
p
p
r
o
x
i
m
a
t
e
t
h
e
J
a
c
o
b
i
a
n
i
f
t
h
e
i
n
p
u
t
J
i
s
e
m
p
t
y
(
u
s
e
i
i
s
s
e
e
m
m
p
p
t
t
y
y
(
(
J
J
)
)
t
o
c
h
e
c
k
)
.
T
o
p
r
o
v
i
d
e
a
n
e
m
p
t
y
i
n
p
u
t
,
u
s
e
[
[
]
]
i
n
t
h
e
c
a
l
l
.
How many steps are required with the analytical Jacobian compared
to the numerical Jacobian? Which algorithm takes longer?
16
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
A
s
s
i
g
m
e
n
t
2
:
R
a
d
i
o
a
c
t
i
v
e
d
e
c
a
y
A decaying radioactive element changes its concentration
according to the ODE
where
k
is the decay constant
The analytical solution to this problem reads
17
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
A
s
s
i
g
n
m
e
n
t
2
(
c
o
n
t
i
n
u
e
d
)
1.
Plot the behavior of the radioactive decay problem as a
vector field in the y vs t plane
Find online the function vector_field.m. It plots the solutions and
derivatives of first order IVPs for different initial conditions
Plot the vector field for different initial values between 0 and 1, time
between 0 and 2 and k = 1. Can you see what a solver has to do?
Is it possible to switch from one trajectory in the vector field to
another? What follows for the uniqueness of the solutions?
2.
The forward Euler method reads for this problem
First define a function defining the successor of y
n
with a header like
function
ynp = Forward_stepper(k, yn, h)
ynp = Forward_stepper(k, yn, h)
Use the conditions y0 = 1, k = 1 and h = 0.1 to solve the radioactive
decay problem from t0 = 0 to tEnd = 10 by defining an integrating
function of the form
function [T,Y] =
function [T,Y] =
stepper_integrate(@stepper,k,t0,tEnd,y0,h)
stepper_integrate(@stepper,k,t0,tEnd,y0,h)
18
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
A
s
s
i
g
n
m
e
n
t
2
(
c
o
n
t
i
n
u
e
d
)
19
3.
The backward Euler method uses the following step formula
(5)
a.
Rearrange this equation so that you can define a function defining the successor of yn
with a header like function ynp = Backward_stepper(k, yn, h)
b.
Use the prior defined integration function to solve the problem with the backward
method.
4.
Plot the obtained solutions comparing them to the analytical solution.
Use the subplot method to produce two subplots in a single figure. What
happens if you increase the step size h? What happens if you increase
the step size above 2?
Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers
Consider the fundamentals of Jacobian matrices, zero of nonlinear equation systems, function linearization, and the multidimensional Newton algorithm in the context of explicit ODE solvers for chemical engineering applications. Dive into the theoretical concepts and practical implementations presented by Sarah Duclos Ivetich from ETH Zurich, focusing on finding solutions to complex equations and equations systems using numerical techniques.
Download Presentation
Please find below an Image/Link to download the presentation.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.
E N D
Presentation Transcript
Jacobian Matrices Catch Up Sarah Duclos Ivetich ETH Zurich, Institut f r Chemie- und Bioingenieurwissenschaften ETH H nggerberg / HCI F106 Z rich E-Mail: sarah.duclos@chem.ethz.ch https://shihlab.ethz.ch/education/Snm/numerics.html Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 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 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 Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 2
Function linearization ( ( ) ) , , f x y f x y ( ) f x 1 = = 0 2 Let us again consider Taylor expansion f x f y + + = ( , ) f x y ( , f x y ) ( ) ( ) 0 1 x x 1 y y 1 1 0 0 0 0 ( , ) x y ( , ) x y 0 0 0 0 f x f y + + = ( , ) f x y ( , f x y ) ( ) ( ) 0 2 x x 2 y y 2 2 0 0 0 0 ( , ) x y ( , ) x y 0 0 0 0 In matrix form, this reads / / f x ( , ( , f x y ) ) / / x y x y f x y f x f f y y 0 1 0 0 1 1 = ( ) 0 2 0 0 2 2 , x y 0 0 Which is equivalent to ? ?? ???= ?(?(?)) Newton method ?(?+1)= ?(?)+ ??(?) Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 3
The Jacobian matrix Consider a continuous, differentiable function ?: ? ? ( ( ) ) , , , , , , , , f x x x f x x x x x 1 1 2 3 n f x f x f x f x f x f x 2 1 2 3 n 1 1 1 ( ) , , , , f x x x x 1 2 n 1 2 3 n n 2 2 2 f x = J = i J , i k 1 1 n k f x f x f x n n n 1 2 n Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 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 ? ?????= ?(?(?)) ?(?+1)= ?(?)+ ??(?) Iterate 2 through 4 until some stopping criteria are fulfilled Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 5
Ordinary Differential Equations (ODEs) d ( ) d y t = ( , ) f t y t Sarah Duclos Ivetich ETH Zurich, Institut f r Chemie- und Bioingenieurwissenschaften ETH H nggerberg / HCI F106 Z rich E-Mail: sarah.duclos@chem.ethz.ch https://shihlab.ethz.ch/education/Snm/numerics.html Sarah Duclos Ivetich / Numerical Methods for Chemical Engineers / Explicit ODE Solvers 6
Problem Definition We are going to look at initial value problems of explicit ODE systems, which means that they can be cast in the following form d ( ) d ( ) y t y t = ( , ) f t y t = y 0 0 Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 7
Example: A 1-D First Order ODE Consider the following initial value problem ( ) 2 y t dt y = dy t = 1 4 ( ) = + + 2 t ( ) y t 3e 2 1 t (0) 1 7 This ODE and its solution follow a general form (Johann Bernoulli, Basel, 17th 18thcentury) 6 5 d ( ) d y t = ( ) r t ( ) ( ) p t y t 4 y t 3 0.1 ( )d p s ( )d p s s s = + ( ) y t e ( )e r t d k t 2 1 0 0.2 0.3 0.4 0.5 t 0.6 0.7 0.8 0.9 1 Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 8
Some Nomenclature On the next slides, the following notation will be used Value of at point Value of the derivative of at point Integration step y y h y t n n y t n n Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 9
Numerical Solution of a 1-D First Order ODE ( ) dt dy t = y t 2 1 4 ( ) = + + 2 t ( ) y t 3e 2 1 t = (0) 1 y 7 Since we know the first derivative, Taylor expansion suggests itself as a first approach 6 5 ( ) d d y t ( ) 2 + + ( ) y y t t O t t 4 x + + + 1 1 1 n n n n n n t 3 n + ( ) O h ( ) O h Exact Solution n = + = h f t y + + 2 2 ( , n ) y hy y 2 n n n Euler Method 1 This is known as the forward / explicit Euler method, named after Leonhard Euler (Basel, 18thcentury) 0 0.1 0.2 0.3 0.4 0.5 t 0.6 0.7 0.8 0.9 1 Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 10
Definition of an implicit algorithm In an implicit algorithm, the function value at the next step appears on the right hand side of the step equation: = F ( , , , ) y t y y + + + 1 1 1 n n n n A simple example is the backward / implicit Euler method: ( , n n n n y y hf t y + + = + ) + 1 1 1 A system of equations has to be solved at every iteration step; Why even use implicit algorithms? Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 11
How does Matlab do it? Non-Stiff Solvers ode45: Most general solver, best to try first; Uses an explicit Runge-Kutta pair (Dormand-Prince, p = 4 and 5) ode23: More efficient than ode45 at crude tolerances, better with moderate stiffness; Uses an explicit Runge- Kutta pair (Bogacki-Shampine, p = 2 and 3). ode113: Better at stringent tolerances and when the ODE function file is very expensive to evaluate; Uses a variable order Adams-Bashforth-Moulton method Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 12
How does Matlab do it? Stiff Solvers ode15s: Use this if ode45 fails or is very inefficient and you suspect that the problem is stiff; Uses multi-step numerical differentiation formulas (approximates the derivative in the next point by extrapolation from the previous points) ode23s may be more efficient than ode15s at crude tolerances and for some problems (especially if there are largely different time-scales / dynamics involved); Uses a modified Rosenbrock formula of order 2 (one-step method) Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 13
Matlab Syntax Hints All ODE-solvers use the syntax [T, Y] = ode45(ode_fun, tspan, y0, ); ode_fun is a function handle taking two inputs, a scalar t (current time point) and a vector y (current function values), and returning as output a vector of the derivatives dy/dt in these points tspan is either a two component vector [tstart, tend] or a vector of time points where the solution is needed [tstart, t1, t2, ...] y0 are the values of y at tstart Use parametrizing functions to pass more arguments to ode_fun if needed f_new = @(t,y)ode_fun(t, y, k, p); This can be done directly in the solver call [T, Y] = ode45(@(t,y)ode_fun(t, y, k, p), tspan, y0); Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 14
Assignment 1 1. Write down the analytical Jacobian matrix for the steady state CSTR; 2. Implement the basic Newton method as shown. 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. Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 15
Assignment 1 (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? Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 16
Assigment 2: Radioactive decay A decaying radioactive element changes its concentration according to the ODE d d y t= ky where k is the decay constant The analytical solution to this problem reads ( ) ( ) y t y kt = 0exp Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 17
Assignment 2 (continued) 1. Plot the behavior of the radioactive decay problem as a vector field in the y vs t plane Find online the function vector_field.m. It plots the solutions and derivatives of first order IVPs for different initial conditions Plot the vector field for different initial values between 0 and 1, time between 0 and 2 and k = 1. Can you see what a solver has to do? Is it possible to switch from one trajectory in the vector field to another? What follows for the uniqueness of the solutions? 2. The forward Euler method reads for this problem ( 1 , + n n n y y h f t y ) += = h k y y n n n First define a function defining the successor of yn with a header like function ynp = Forward_stepper(k, yn, h) Use the conditions y0 = 1, k = 1 and h = 0.1 to solve the radioactive decay problem from t0 = 0 to tEnd = 10 by defining an integrating function of the form function [T,Y] = stepper_integrate(@stepper,k,t0,tEnd,y0,h) Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 18
Assignment 2 (continued) The backward Euler method uses the following step formula 3. (5) a. Rearrange this equation so that you can define a function defining the successor of yn with a header like function ynp = Backward_stepper(k, yn, h) b. Use the prior defined integration function to solve the problem with the backward method. Plot the obtained solutions comparing them to the analytical solution. Use the subplot method to produce two subplots in a single figure. What happens if you increase the step size h? What happens if you increase the step size above 2? 4. Sarah Duclos Ivetich/ Numerical Methods for Chemical Engineers / Explicit ODE Solvers 19
