1 / 25

Ordinary Differential Equations (ODEs)

Ordinary Differential Equations (ODEs). Daniel Baur ETH Zurich, Institut für Chemie- und Bioingenieurwissenschaften ETH Hönggerberg / HCI F128 – Zürich E-Mail: daniel.baur@chem.ethz.ch http://www.morbidelli-group.ethz.ch/education/index . Problem Definition.

shona
Download Presentation

Ordinary Differential Equations (ODEs)

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

More Related