Page 222 - Numerical methods for chemical engineering

208 4 Initial value problems

the ODE system ˙x = f (t, x; Θ), we must at a minimum supply a routine that returns
f (t, x; Θ) with the syntax
function f = calc f(t, x, P1, P2, . . .);
P1, P2, . . . are optional fixed parameters. The IVP is solved by

[t traj, x traj] = ode45(@calc f, t span, x 0, Options, P1, P2, . . .);
t span contains the start and end times of the simulation, or optionally a list of discrete times
at which the states are needed. x 0 is the initial state and Options is a structure, managed
by odeset, that allows the user to modify the solver behavior (use [ ] to accept the default
options). P1, P2, . . . are optional fixed parameters passed to the function routine, here
calc f.
T
Implicit methods such as ode15s require the Jacobian matrix J = ∂f /∂x . If only a
routine such as calc f is supplied, ode15s estimates the Jacobian by finite differences.
For large, sparse systems, the work associated with Jacobian estimation can be reduced
significantly by supplying a matrix with the same sparsity pattern as the Jacobian, through
the “JPattern” field of Options. Even better, we can supply a routine that computes the
Jacobian,
function Jac = calc Jac(t, x, P1, P2, . . .);
and set the “Jacobian” field of Options to its name, here ‘calc Jac’.
For a DAE system M ˙ x = f (t, x; Θ), ode15s may be used if the system is of index one,
supplying either a constant mass matrix or a function,
function M = calc Mass(t, x);
and informing ode15s of its identity through the “Mass” field of Options, along with the
appropriate specifications of “MStateDependence” and “MassSingular.” More general
DAE systems of the fully implicit form F(t, x, ˙ x) = 0 can be solved by ode15i. The low-
order implicit methods ode23s and ode23tb can be useful for very stiff systems, such as
discretized PDEs (discussed later in Chapter 6).
odeplot plots the state trajectory, and odeprint prints it to the screen. Related functions
for 2-D and 3-D plots are odephas2 and odephas3.

Problems

4.A.1. Consider the following ODE system with a steady state at (1, 2),
2
˙ x 1 = f 1 (x 1 , x 2 ) =−2(x 1 − 1)(x 2 − 2) + (x 1 − 1)(x 2 − 2) − (x 2 − 2) − 3(x 1 − 1)
2
˙ x 2 = f 2 (x 1 , x 2 ) =−(x 1 − 1) + (x 1 − 1) (x 2 − 2) (4.247)
Derive the linearized ODE system that describes the response of the system to small per-
turbations around the steady state. Is the steady state stable? Will the steady state exhibit
an oscillatory response to small perturbations? Confirm your results by writing a program
to simulate the response to a random, small perturbation.