ME451 Kinematics and Dynamics of Machine Systems

1. ME451 Kinematics and Dynamics of Machine Systems Numerical Integration. Stiffness: Implicit Methods. October 30, 2013 Radu Serban University of Wisconsin-Madison

2. Before we get started… • Last Time: • Recovering constraint reaction forces • Started discussion on numerical integration • Today • Stiff differential equations • Implicit numerical integration formulas • Assignments: • Homework 9 – 6.3.3, 6.4.1 – due November 4 (12:00pm) • Matlab 7 – due tonight, October 30, Learn@UW (11:59pm) • Project 1 – due Wednesday, November 6, Learn@UW (11:59pm) • Your simEngine2D can now perform complete Kinematic analysis! • Test the provided visualization function (available on the SBEL course webpage) • Miscellaneous • Draft proposals for the Final Project due on Friday, November 1 • Midterm 2 – Wednesday, November 6, 12:00pm in ME 1143 • Review session – Monday, November 4, 6:30pm in ME 1143

3. Numerical Integration

4. Basic Concept • IVP • In general, all we can hope for is approximating the solution at a sequence of discrete points in time • Uniform grid (constant step integration) • Adaptive grid (variable step integration) • The numerical solution is then the sequence of approximations • Basic idea: somehow turn the differential problem into an algebraic problem (approximate the derivatives)

5. Simplest method: Forward Euler • Starting from the IVP • Use the simplest approximation to the derivative • Rewrite the above asand use ODE to obtain Forward Euler Method with constant step-size

6. FE: Geometrical Interpretation • IVP • Forward Euler integration formula

7. FE: Example • Solve the IVPusing Forward Euler (FE) with a step-size • Compare to the exact solution

8. Forward Euler: Effect of Step-Size % IVP (RHS + IC) f = @(t,y) -0.1*y + sin(t); y0 = 0; tend = 50; % Analytical solution y_an = @(t) (0.1*sin(t) - cos(t) + exp(-0.1*t)) / (1+0.1^2); % Loop over the various step-size values and plot errors colors = [[0, 0.4, 0]; [1, 0.5, 0]; [0.6, 0, 0]]; Figure, hold on, box on h = [0.001 0.01 0.1]; for ih = 1:length(h) tspan = 0:h(ih):tend; y = zeros(size(tspan)); err = zeros(size(tspan)); y(1) = y0; err(1) = 0; for i = 2:length(tspan) y(i) = y(i-1) + h(ih) * f(tspan(i-1), y(i-1)); err(i) = y(i) - y_an(tspan(i)); end plot(tspan, err, 'color', colors(ih,:)); end legend('h = 0.001', 'h = 0.01', 'h = 0.1'); FE errors for different values of the step-size

9. FE: Effect of Step-Size % IVP (RHS + IC) f = @(t,y) -0.1*y + sin(t); y0 = 0; tend = 50; % Loop over the various step-size values and plot errors colors = [[0, 0.4, 0]; [1, 0.5, 0]; [0.6, 0, 0]]; Figure, hold on, box on h = [0.1 1.0 5.0]; for ih = 1:length(h) tspan = 0:h(ih):tend; y = zeros(size(tspan)); y(1) = y0; for i = 2:length(tspan) y(i) = y(i-1) + h(ih) * f(tspan(i-1), y(i-1)); end plot(tspan,y, 'color', colors(ih,:)) end legend('h = 0.1', 'h = 1', 'h = 5'); FE accuracy at different values of the step-size

10. Stiff Differential Equations

11. A Simple IVP Example • Consider the IVPwhose analytical solution is • 5 integration steps with Forward Euler formula, , , • Compare the global error (difference between analytical and numerical solution) Error for 0 0.36787944117144 0.13533528323661 0.04978706836786 0.01831563888873 0.00673794699909 Error for 0 2.04978706836786 -3.99752124782333 8.00012340980409 -15.99999385578765 32.00000030590232 Error for 0 0.01873075307798 0.03032004603564 0.03681163609403 0.03972896411722 0.04019944117144

12. A Simple IVP Example Analytical Solution Forward Euler

13. Stiff Differential Equations • Problems for which explicit integration methods (such as Forward Euler) do not work well • Other explicit formulas: Runge-Kutta (RK4), DOPRI5, Adams-Bashforth, etc. • Stiff problems require a different class of integration methods: implicit formulas • The simplest implicit integration formula: Backward Euler (BE)

14. BE: Geometrical Interpretation • IVP • Forward Euler integration formula • Backward Euler integration formula

15. A Simple IVP Example Analytical Solution Backward Euler Forward Euler

16. Forward Euler vs. Backward Euler

17. Backward-Difference Formulas (BDF) • IVP • Also known as Gear method (DIFSUB – 1971) • Family of implicit linear multi-step formulas • BDF of 1st order: • BDF of 2nd order: • BDF of 3rd order: • BDF of 4th order: • BDF of 5th order: C.W. (Bill) Gear b. 1935

18. Curtiss & Hirschfelder Example Forward Euler Backward Euler C.F. Curtiss and J.O. Hirschfelder – “Integration of Stiff Equations”Proc. Nat. Acad. Sci, USA (1952)

19. Stability and Accuracy of Numerical Integrators

20. Two Key Properties of a Numerical Integrator Two properties are relevant when considering a numerical integrator for finding an approximation of the solution of an IVP • Stability • Accuracy

21. Stability of a Numerical Integrator • The problem: • How big can the integration step-size be without the numerical solution blowing up? • Tough question, answered in a Numerical Analysis class • Different integration formulas, have different stability regions • We would like to use an integration formula with large stability region: • Example: Backward Euler, BDF methods, Newmark, etc. • Why not always use these methods with large stability region? • There is no free lunch: these methods are implicit methods that require the solution of an algebraic problem at each step

22. Accuracy of a Numerical Integrator • The problem: • How accurate is the formula that we are using? • If we decrease , how will the accuracy of the numerical solution improve? • Tough question, answered in a Numerical Analysis class • Examples: • Forward and Backward Euler: accuracy • RK45: accuracy • Why not always use methods with high accuracy order? • There is no free lunch: these methods usually have very small stability regions • Therefore, they are limited to using very small values of

23. Implicit Integration of Stiff ODEs

24. Implicit Integration and Nonlinear Systems[Examples] • Backward Euler (BE) formula: • Apply one step of BE for the following IVPs ()

25. Implicit Integration and Nonlinear Systems[Example 1] • Backward Euler (BE) formula: • Apply one step of BE for the following IVP () • This example shows that discretization with an implicit integration formula requires the solution of an algebraic equation

26. Implicit Integration and Nonlinear Systems[Example 2] • Backward Euler (BE) formula: • Apply one step of BE for the following IVP () • This example shows that the solution of the resulting algebraic equation raises additional problems (selecting one of multiple possible solutions)

27. Implicit Integration and Nonlinear Systems[Example 3] • Backward Euler (BE) formula: • Apply one step of BE for the following IVP () • This example shows that, in general, the resulting nonlinear algebraic problem can only be solved numerically.

28. Implicit Integration and Nonlinear Systems[Examples] • Backward Euler (BE) formula: • Apply one step of BE for the following IVPs () • Upon discretization with an implicit integration formula (such as BE), the approximation to the solution of the DE at the new time is obtained by solving an algebraic equation • In general, the resulting (system of) nonlinear equation(s) can only be solved numerically (using for example the Newton-Raphson method)

29. Implicit Integration: Conclusions • Stiff problems require the use of implicit integration methods • Because they have very good stability, implicit integration methods allow for step-sizes that could be orders of magnitude larger than those needed if using explicit integration • However, for most real-life IVPs, discretization using an implicit integration formula leads to another nasty problem: • To find the solution at the new time, we must solve a nonlinear algebraic problem • This brings back into the picture the Newton-Raphson method (and its variants) • We have to deal with providing a good starting point (initial guess), computing the matrix of partial derivatives, etc.

30. Solution of IVPs in Matlab

31. General Call to a Matlab ODE Solver • IVP: • Example of using Matlab’sode45: • Use odesetto create an ODE options structure to change default values for various solver parameters [t,y] = ode45(odefun,tspan,y0,options) function dydt = odefun(t,y) y0 [t0 tend] options = odeset(…)

32. Example of IVP Solution with Matlab Convert to 1st order ODE >> [t,y]=ode45(‘rhs’,[0 20],[2;0]); >> plot(t, y); Convert tovector form function yd = rhs(t, y) yd = [y(2); (1-y(1)^2)*y(2)-y(1)];

33. ODE Solvers in MATLAB Use ode45 for non-stiff problems and ode15s for stiff problems

34. A Stiff Problem ---- ode45 solver statistics ---- 7557 successful steps 504 failed attempts 48367 function evaluations Elapsed time is 0.564820 seconds. ---- ode15s solver statistics ---- 139 successful steps 3 failed attempts 288 function evaluations 1 partial derivatives 27 LU decompositions 284 solutions of linear systems Elapsed time is 0.188120 seconds. % specify RHS function rhs = @(t,y) [-1,-1;1,-5000]*y; % specify accuracy and turn on stats option options = odeset('RelTol',1e-6,'Stats','on'); % call the ode45 solver fprintf('\n---- ode45 solver statistics ----\n') tic, [~,~] = ode45(rhs, [0, 5] , [1; 1], options); toc % call the ode15s solver fprintf('\n---- ode15s solver statistics ----\n') tic, [~,~] = ode15s(rhs, [0, 5] , [1; 1], options); toc • On a stiff problem, • an explicit solver is forced to take small steps to ensure stability • an implicit solver can take much larger steps, which reduces overall time to solution; but it must solve algebraic equations at each time step