1 / 37

ME451 Kinematics and Dynamics of Machine Systems

ME451 Kinematics and Dynamics of Machine Systems. Dealing with Differential Equations ~ A Numerical Approach ~ March 12, 2009. © Dan Negrut, 2009 ME451, UW-Madison. Before we get started…. Last Time Exam Before that: talked about singularities of mechanism

emily
Download Presentation

ME451 Kinematics and Dynamics of Machine Systems

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. ME451 Kinematics and Dynamics of Machine Systems Dealing with Differential Equations~ A Numerical Approach ~ March 12, 2009 © Dan Negrut, 2009ME451, UW-Madison

  2. Before we get started… • Last Time • Exam • Before that: talked about singularities of mechanism • Remember that as long as the constraint Jacobian is noinsingular, you are in great shape • Today • HW assigned (due on March 26): • Problems stated on two separate slides of today’s lecture • Basically require you to solve an IVP using MATLAB • ADAMS component has been uploaded • We’ll discuss about • Numerical Methods • Initial Value Problems and their numerical solution • Material not included in the book 2

  3. Numerical Method(also called Numerical Algorithm, or simply Algorithm) • Represents a recipe, a succession of steps that one takes to find the solution of a problem that does not admit an analytical solution • Analytical solution: sometimes called “closed form”, or “exact” solution • The solution obtained with the numerical method is called “numerical solution” • Examples: • Evaluate the integral • Solve the equation • Solve the IVP that comes form the EOM for a simple pendulum • Many, many others (actually very seldom can you find the exact solution of a problem…) 3

  4. Where are Numerical Methods Used? • Thermo and Heat Transfer: compute the flow of heat (diffusion) through a media (Finite Element Methods) • Fluid Dynamics: solve the Navier-Stokes equations (Finite Element/Finite Volume Methods) • Mechanics of Materials: solve the partial differential equations to find stress/strain distribution (Finite Element Methods) • Multibody Dynamics: solve differential-algebraic equations (DAEs) that result from Newton’s second law 4

  5. Numerical Methods in ME451 • In regards to ME451, one would use numerical method to solve the dynamics problem (the resulting set of differential equations that capture Newton’s second law) • The particular class of numerical methods used to solve differential equations is typically called “numerical integrators”, or “integration formulas” • A numerical integrator generates a numerical solution at discrete time points (also called grid points, station points, nodes) • This is in fact just like in Kinematics, where the solution is computed on a time grid • Different numerical integrators generate different solutions, but the solutions are typically very close together, and [hopefully] closed to the actual solution of our problem • In 99% of the cases, the use of numerical integrators is the only alternative for solving complicated systems described by non-linear differential equations 5

  6. Numerical Integration~Basic Concepts~ • Initial Value Problem: (IVP) • So what’s the problem? • You are looking for a function y(t) that depends on time (changes in time), whose time derivative is equal to a function f(t,y) that is given to you (see IVP above) • In other words, I give you the derivative of a function, can you tell me what the function is? • In ME451, the best you can hope for is to find an approximation of the unknown function y(t) at a sequence of discrete points (as many of them as you wish) • The numerical algorithm produces an approximation of the value of the unknown function y(t) at the each grid point. That is, the numerical algorithm produces y(t1), y(t2), y(t3), etc. 6

  7. Relation to ME451 • When carrying out Dynamics Analysis, what you can compute is the acceleration of each part in the model • Acceleration represents the second time derivative of your coordinates • Somewhat oversimplifying the problem to make the point across, in ME451 you get the second time derivate • This represents a second order differential equation since it has two time derivatives taken on the position q • Problem is reduced to a set of first order differential equations by introducing a helper variable v (the velocity): • With this, the original second order differential problem becomes: 7

  8. The difference between an ordinary differential equation (ODE) and an initial value problem (IVP) 8

  9. ODE vs. IVP • Difference between ODE and IVP - Often a source of confusion • Ordinary Differential Equation (ODE) • Typically, has an infinite number of solutions • Initial Value Problem (IVP) • Is an ODE plus an initial condition (IC): • The solution assumes at time T=0 a certain prescribed value • The solution for the IVP’s that we’ll deal with is UNIQUE 9

  10. ODE: Infinite Number of Solutions • ODE Problem: • Initial Condition: y0 =[-1000:50:1000] 10

  11. ODE vs. IVP (Contd.) • Example: consider a simple first order ODE • It has an infinite number of solutions: • However, if you specify an Initial Condition (IC) at time t=0, for instance, x(0) =2.5, then there is a UNIQUE solution that satisfies both the ODE and the imposed IC: • Remember: 1. IVP = ODE + IC 2. IVP has a UNIQUE solution (unlike on ODE) 11

  12. How do you go about solving an Initial Value Problem? A look at: Euler’s Method Predictor-Corrector Method Runge-Kutta Method 12

  13. Find solution of this Initial Value Problem: Numerical Integration:Euler’s Method Euler’s Method: • t is the step size • The idea: at each grid point tk, turn the differential problem into an algebraic problem by approximating the value of the time derivative: 13

  14. Example:- Integrate 5 steps using Euler’s Method- Compare to exact solution f(t,y) = -10y (note no explicit dependency on time t for f) k=0 y0=1.0 k=1 y1=y0+f(t0,y0)t = 1.0000 + (-10*1.000)*0.01 = 0.9000 k=2 y2=y1+f(t1,y1)t = 0.9000 + (-10*.9000)*0.01 = 0.8100 k=3 y3=y2+f(t2,y2)t = 0.8100 + (-10*.8100)*0.01 = 0.7290 k=4 y4=y3+f(t3,y3)t = 0.7290 + (-10*.7290)*0.01 = 0.6561 k=5 y5=y4+f(t4,y4)t = 0.6561 + (-10*.6561)*0.01 = 0.5905 Exact solution: y(t) = e-10t Solution: y(0) =1.0000 y(0.01)=0.9048 y(0.02)=0.8187 y(0.03)=0.7408 y(0.04)=0.6703 y(0.05)=0.6065 14

  15. Example- Integrate in Matlab for 1 second using Euler’s Method- Compare to exact solution • Matlab y0=1; % Euler’s Method dt=.01; t=0:.01:1.; yh=y0; for i=2:length(t) f=-10*yh(i-1); yh(i)=yh(i-1)+f*dt; end % Exact solution y=y0*exp(-10*t); % Plot and compare plot(t,y,'b-',t,yh,'ro'); 15

  16. Euler Method: ~ Effect of Step Size ~ • Solve using step sizes t=0.1, 1 and 5 sec dt=5 sec • Matlab dt=1 sec dt=0.1 sec % Exact solution t=0:.01:50; y=.99*exp(-t/10)+.995*sin(t-1.47); % Numerical solution dt=.1; th=0:dt:50; yh=0; for i=2:length(th) f=-yh(i-1)/10+sin(th(i)); yh(i)=yh(i-1)+f*dt; end plot(t,y,'b-',th,yh,'ro-'); Conclusion: If you use large step-sizes  t, the quality of the solution is really bad (you can’t be too aggressive) 16

  17. Predictor-Corrector Methods • Instead of computing f(y,t) at one point, you can average over multiple points (two in this case) • To implement, must predict yk+1 using forward Euler Method The predictor The corrector 17

  18. Runge-Kutta Methods • 4th Order Runge-Kutta • Uses slope at 4 pts within each time step • Per-step error is proportional to t5 • Euler’s Method [introduced a couple of slides ago…] • Simple first order • Per-step error is proportional to t2 18

  19. MATLAB Support for solving IVP 19

  20. Ordinary Differential Equations(Initial Value Problem) • An ODE + initial value: • Use ode45 for non-stiff IVPs and ode23t for stiff IVPs (concept of “stiffness” discussed later in presentation) [t,y] = ode45(odefun,tspan,y0,options) function dydt = odefun(t,y) Initialvlue [initialtime finaltime] • Use odeset to define options parameter 20

  21. IVP Example: function dydt = myfunc(t,y) dydt=zeros(2,1); dydt(1)=y(2); dydt(2)=(1-y(1)^2)*y(2)-y(1); • [t,y]=ode45('myfunc',[0 20],[2;0]) Note: Help on odesetto set options for more accuracy and other useful utilities like drawing results during solving. 21

  22. Homework Problem 1 (out of 2) • Use MATLAB’s ODE45 in order to solve the IVP below on the interval [0,20] seconds: 22

  23. Homework Problem 2 (out of 2) • Consider the following IVP: • In MATLAB, over the interval t 2 [0,100]: • Implement the Backward Euler (BE) formula to find the solution of this IVP • Work with a step-size  t that ensures good quality results • Implement the Runge-Kutta (RK) method discussed to solve the IVP • Work with a step-size  t that ensures good quality results • Use MATLAB’s ODE15s numerical integrator to solve the IVP • Find the analytical solution of the above IVP and compare it with the numerical solutions obtained above (that is, plot all solutions together) • Of BE, RK, and ODE15s, which seems to finish first when computing the solution? Does this come in line with expectations and why? 23

  24. The concept of stiff differential equations, and how to solve the corresponding IVP 24

  25. Example: IVP 0 0.01873075307798 0.03032004603564 0.03681163609403 0.03972896411722 0.04019944117144 - Integrate 5 steps using forward Euler formula: t=0.002, t=0.01, t=0.03 - Compare the errors between numerical and analytical solutions t=0.002: Error t=0.01:Error t=0.03:Error 0 0.36787944117144 0.13533528323661 0.04978706836786 0.01831563888873 0.00673794699909 0 2.04978706836786 -3.99752124782333 8.00012340980409 -15.99999385578765 32.00000030590232 25

  26. Example: Analytical Solution Forward Euler (t=0.03) 26

  27. Concept of stiff IVP’s • IVP’s for which forward Euler doesn’t work well (see example) • In general, the entire class of so called explicit formulas doesn’t work • Forward Euler, Runge-Kuta (RK23, RK45), DOPRI5, etc. • Stiff IVP’s require a different class of integration formulas • Implicit formulas • Example: backward Euler 27

  28. Forward Euler • Backward Euler Explicit vs. Implicit Formulas(look at Euler family) • Initial Value Problem 28

  29. The Nonlinear System Associated with an Implicit Integrator • Implicit integration requires finding the value of yk+1 that verifies the equation • You find it by solving the following nonlinear equation • The nonlinear function  is defined as: 29

  30. Back to Newton-Raphson… • Newton algorithm for nonlinear systems requires: • An initial guess from where solution starts being searched for • An iterative process in which the approximation of the solution is gradually improved: • I’m dropping the subscript “k+1” to keep notation simple. Then, 30

  31. Example: Exact Solution Backward Euler and Analytical Solution Forward Euler (t=0.03) 31

  32. Stiffness: Detecting It. • The IVP we just considered was stiff: • Dealing with stiff IVP’s is tricky • Calls for special integration formulas (we just saw that) • How can I tell if a system is stiff ? • Two ways: • You have “inside” information about the problem you solve (lots of damping in a mechanical system, for instance, leads to stiffness) • You compute and analyze the eigenvalues of the “Jacobian” matrix. If eigenvalues are way out to the left of the complex plane, the system is stiff… 32

  33. Numerical Integration Formula:The STABILITY Attribute • The stability question: • So how big can I choose the integration step-size t and be safe? • Tough question (Numerical Analysis class) • Different integration formulas, have different stability regions • You’d like to use an integration formula with large stability region: • Example: Backward Euler, BDF methods, Newmark, HHT • 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 algebra problem at each step 33

  34. Numerical Integration Formula :The ACCURACY Attribute • The accuracy question: • How accurate is the formula that I’m using? • If I start decreasing t , how will the accuracy of the numerical solution scheme improve? • Tough question (Numerical Analysis class) • Examples: • Forward and Backward Euler: accuracy O(t ) • RK45: accuracy O(t 4) • Why not always use methods with high accuracy order? • There is no free lunch: these methods usually have very small stability • Therefore, you are limited to very small values of t 34

  35. ODE solvers in Matlab 35

  36. BDF Integrators(Used for Stiff IVPs) • 1st order: • 2nd order: • 3rd order: • 4th order: • 5th order: 36

  37. Concluding Remarks • The purpose of a numerical method is to find an approximation of the solution of a mathematical problem when you cannot find an analytical (closed form, exact) solution • Powerful and inexpensive computers revolutionized the use of numerical methods • Simulation of a car crash in minute detail • Formation of galaxies • Folding of a protein • Finding the electron distribution around nuclei in a nanostructure • Numerical methods enable the concept of “simulation-based engineering” • You use computer simulation to understand how your mechanism (mechanical system) behaves, how it can be modified and controlled 37

More Related