Ordinary Differential Equations (ODEs)

Ordinary Differential Equations (ODEs) Daniel BaurETH Zurich, Institut für Chemie- und BioingenieurwissenschaftenETH Hönggerberg / HCI F128 – ZürichE-Mail: daniel.baur@chem.ethz.chhttp://www.morbidelli-group.ethz.ch/education/index Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

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 Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Higher Order ODEs • Higher order explicit ODEs can be cast into the first order form by using the following «trick» • Therefore, a system of ODEs of order n can be reduced to first order by adding n-1 variables and equations Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Autonomous Systems • A system is called autonomous if the independent variable (t) does not explicitly appear in the system • It is always possible to transform a non-autonomous system into an autonomous one Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Example: A 1-D First Order ODE • Consider the following initial value problem • This ODE and its solution follow a general form(Johann Bernoulli, Basel, 17th– 18th century) Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Numerical Solution of a 1-D First Order ODE • Since we know the first derivative, Taylor expansion suggests itself as a first approach • This is known as the Euler method, named after Leonhard Euler (Basel, 18th century) Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Some Nomenclature • On the next slides, the following notation will be used Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Categorization of Algorithms • One step integration algorithms: • Multi-step integration algorithms: • Explicit integration algorithms: • Implicit integration algorithms: Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Accuracy of Integration Algorithms • The local error is always proportional to a known power of the integration step • An algorithm is of orderp if it can correctly integrate a polynomial of order p • An algorithm of order p always has a local error in the order of hp+1 • An algorithm is called convergent if for h  0 the global error decreases (assuming there are no rounding errors) Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Examples • Forward Euler Algorithm • Explicit Runge-Kutta (p = 2) h = 0.01h = 0.1h = 0.2 h = 0.01h = 0.2 • Error at t = 1 (h = 0.01) • Euler Forward: 1.72% • RK (p=2): 0.01% Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Conditioning of the Problem • A system of ODEs is called well conditioned if a small error in the function (e.g. parameters) and/or in the initial conditions generate a small error in the solution (without rounding errors!) • Example of perturbations / errors: • The conditioning of a problem has nothing to do with the stability of an algorithm! Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Example of Bad Conditioning • Consider the ODE • The general solution reads • If the initial condition y(0) = 1 applies, the constant is 0 • On the other hand, if y(0) = 1.0001, the solution reads Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Examples for Conditioning • Consider these initial value problems • They all have the same solution • But if we perturb the initial condition (y(0) = 1.0001) Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Measuring the Conditioning • We can tell a priori if a system is ill conditioned:If all the eigenvalues of the Jacobian of a system of ODEs have a positive real part, then the system is ill conditioned. Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Stability of the Algorithms • The Forward Euler algorithm is stable only ifwhere κ is the amplification factor • Therefore • If fy > 0 (the ODE is ill conditioned)  |κ| > 1 for h > 0 • If fy < 0 (the ODE is well conditioned)  |κ| <1 but only if Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Stability and Precision • When dealing with ODE systems, it is common practice to reduce the problem of stability to the ODEwhere λmax is the largest eigenvalue of the Jacobian of the system. Note that eigenvalues are typically complex. • For example, take a look at the Forward Euler method • The method is stable if • The local error of the method is • With the maximum step size Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Explicit Runge-Kutta Methods • A very commonly used method is RK with p = 4 • Advantages • RK is typically efficient and requires small memory • It is easy to change the integration step • Disadvantages • Efficiency is lower than multi-step methods • It is difficult to evaluate the local error Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Multi-Step Algorithms • Multi-step methods use not only the current function value yn to compute the next value yn+1, but also function values at previous times (yn-1, ...) • One example are the Adams-Bashforth methods(here p = 3) • At the start of the integration, only y0 is known, therefore we have to use some other method to calculate y1 and y2 to get the algorithm started (for example the Euler method) Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Multi-Step Algorithms (Continued) • Advatanges • Typically very efficient • Integration step is easy to change • Error estimation is simple • Order of the method is easily changed • It is easy to compute y(t) outside the grid points • Usually, fewer function evaluations are required • ODEs with a non-identity mass matrix A are easily solved Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Stiff Problems • A problem is considered stiff if e.g. the following apply: • A linear constant coefficient system is stiff if all of its eigenvalues have a negative real part and the stiffness ratio is large • Stability requirements, rather than accuracy considerations, constrain the step length • Some components of the solution show much faster dynamics than others • Explicit solvers require an impracticably small h to insure stability, therefore implicit solvers are preferred • This requires solving large (non-)linear algebraic systems, which will be addressed in a later excercise Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

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 ode45at 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 Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

How does Matlab do it? Stiff Solvers • ode15s: Use this if ode45fails 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) Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

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 neededf_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); Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Exercise • Consider a batch reactor where these reactions take place Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Assignment • Using the following values • k1 = 1; k2 = 10; A0 = 1; B0 = 0 • Find the Jacobian of the system. Which are the eigenvalues of the Jacobian? Is the system well conditioned? • Use ode45to solve system. Find the time where 99% of A has been consumed and use this as the final integration time. Plot the results. • Find online templates for the Forward Euler method and the RK method of order 4. • Add the necessary lines to calculate a step, using the formulas on slides 10 (Euler) and 17 (RK4). Which step size h is needed to have a maximum error of 1% compared to ode45 at the end of the integration? What stepsize was used by ode45? • Now change A0 = 0.9; B0 = 0.2; k2 = 1e5 and try ode45 again. What happened? What can you do to resolve this issue? Plot the resulting concentration profiles. Daniel Baur / Numerical Methods for Chemical Engineers / Explicit ODE Solvers

Ordinary Differential Equations (ODEs)