Computation & Simulation

1. Computation & Simulation Assignment Two Numerical methods for solving ODEs Group 5 Xiang Wang IssraaAlhitiPadraigO'Dwyer Tutor: Dr. Conor Brennan

2. Use of Numerical Solvers • 1-Car companies can improve crash safety of their models • This article describes a virtual crash test dummy simulation in which the crash effects on two software generated android models were analysed. The crash simulation resulted in a series of ODEs which represented dynamic effects on these virtual dummies. By solving the ODEs using numerical methods (Runge–Kutta) the authors were able to compute the generalized coordinates, their time derivatives and the body forces of the two models。 • Bud Fox, Leslie S. Jennings, and Albert Y. Zomaya.1999. Numerical computation of differential-algebraic equations for nonlinear dynamics of multibody android systems in automobile crash simulation. • IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 46, NO. 10, OCTOBER 1999. Available from Inspec(EI Engineering Village) database .http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=790496&isnumber=17164 • [Accessed 18 March 2009].

3. Use of Numerical Solvers • 2- Approximation of artificial satellite trajectories and planetary ephemerides. • ODEs are solved using an integration subroutine which generates the body's coordinates and time derivatives and hence produces the planetary ephemerides and satellite trajectories for a time interval. Computer simulation and graphical output indicates the satellite and planetary positions. • Fox, B.; Jennings, L.S.; Zomaya, A.Y.2000.Numerical computation of differential-algebraic equations for the approximation of artificial satellite trajectories and planetary ephemerides. • Transactions of the ASME. Journal of Applied Mechanics, v 67, n 3, 574-80, Sept. 2000 • Available from Inspec(EI Engineering Village) database . • http://www.engineeringvillage2.org/controller/servlet/Controller?SEARCHID=da37721201089cf30M4886prod2data2&CID=quickSearchAbstractFormat&DOCINDEX=6&database=3&format=quickSearchAbstractFormat • [Accessed 18 March 2009].

4. Use of Numerical Solvers • Animation of a body • Modelling each leg segment as a separate pendulum, and calculating its movements continuously. What you see here are not cartoon animations, each move is calculated from first principles, using differential equations. • http://palaeo.gly.bris.ac.uk/dinosaur/Allostick.gif • http://palaeo.gly.bris.ac.uk/dinosaur/animation.html • [Accessed 18 March 2009].

5. David O’Meara, CEO, Havoc:“The quality and quantity of Irish computer science graduates are too low…” Quality? Not according to Peter Brabazon, Dir., Discover Science & Engineering “One of the best in EU for percentage of 3rd level graduates in science, maths & computing” [Keenan,B. 2008. 'INVISIBLE HAND' IS GUIDING STUDENTS IN RIGHT DIRECTION. The Sunday Independent (Ireland) August 24, 2008 O’Meara acknowledges there are plenty of people with ideas but that graduates are more attracted to multinational companies. Why? Something for CEOs of smaller companies to address in terms of incentives for graduates? Possible need to develop business and entrepreneurial skills in graduates Perhaps O’Meara and his cohorts could provide more mentoring and advice rather than criticism of maths and science teaching!

6. David O’Meara, CEO, Havoc:“The quality and quantity of Irish computer science graduates are too low…” • Quantity? • 3rd level computer course intake not up to pre-2002 levels… • Why? • Accepted wisdom: Teaching standards in maths and science and 2nd level facilities are so poor…etc, etc. • But unlikely that standards and facilities have deteriorated so much so quickly! So why the fall off? • Students voting with their feet based on job prospects? Business services predicted as the growth area. [Keenan,B. 2008. 'INVISIBLE HAND' IS GUIDING STUDENTS IN RIGHT DIRECTION. The Sunday Independent (Ireland) August 24, 2008] • Even ESRI Medium Term Review predicts steady fall in manufacturing jobs from 2010 with hi-tech sector shedding 2.5% of jobs per year http://www.esri.ie/UserFiles/publications/20080515155545/MTR11.pdf (p.65-68) • Chicken and egg situation?

7. Projectile Simulation – Matlab Code Objectives Comparison of analytical and numerical methods Adding realism – air friction, bouncing etc. Exploiting Matlab – iterative methods, fast calculations, animation Issues Analytical method => position can be determined at any point in time. BUT computer can only sample time at discrete intervals. Numerical methods. Which one? Euler Predictor-Corrector Runge-Kutta Recursion Accuracy v Speed Code structure Dealing with time variable events – hitting the ground or other object

8. Projectile Simulation – Matlab Code • Code Structure • Setup starting parameters – firing velocity, angle etc. • Step size : Analytical and Numerical • Matrix arrays to store multiple relative values (e.g. x and y co-ords, time value) • WHILE loop • Analytical Method: simply uses equations of motion to calculate the exact position of the projectile at each point in time. x(ct,1)=v_x0*t(ct,1); y(ct,1)=v_y0*t(ct,1)-0.5*g*power(t(ct,1),2); • Numerical Method: successively estimates the position based on the previous estimate and the application of a specific numerical method (Euler in this case) t(ct,1)=t_start+(ct-1)*delta_t; %Calculate velocities in x and y dir numerically if (t(ct,1)~=t_start) v_x(ct,1)=v_x(ct-1,1)+0*delta_t; v_y(ct,1)=v_y(ct-1,1)-g*delta_t; else v_x(ct,1)=v_x0; v_y(ct,1)=v_y0; end %Calculate x and y positions based on velocities x(ct+1,1)=x_0+x(ct,1)+v_x(ct,1)*delta_t; y(ct+1,1)=y(ct,1)+v_y(ct,1)*delta_t; • Loop break out flag – when a certain event occurs if (y(ct+1,1)<threshold)

9. Projectile Simulation – Matlab Code • Plot results • Within loop ensures each iteration plots its result • Repeat loop for different situations • Projectile without air friction applied (analytically and numerically) • Projectile with air friction applied (numerically with different air friction constants)

10. Projectile Simulation – Matlab Code Plot:

11. Caparison of moving ball with/without air friction

12. Challenge One • Runge Kutta is used instead of Euler series • Advantage • Increasing the speed of matlab code • Smaller error when delta time getting bigger and bigger • Euler series : delta time = 0.2s • Elapsed time is 2.150858s • Runge Kutta : delta time = 0.4s • Elapsed time is 1.031713s

13. Runge Kutta vs. Euler Series

14. Challenge Two--hitting

15. Balls Hitting • Two balls being thrown at the same time • Velocities are fixed for both • the throwing angle is fix for upper one • The initial throwing angle is 90 degree for bottom one, it will reduce the angle if it cannot hit the upper ball. • Hitting condition • Time, x distance and y distance should be the same

16. Balls Hitting • Tips in matlab code • Two flags are used in our code b &c • Flags b : speed up of finding the hitting point • If both objects did not collide and either one hits the ground – b=1; • if the objects do collide – b=m+n • If velocity y is great than a number i.e. 1m/s m=0, n=0 • If velocity y is less than 1m/s m=1, n=1 • Flags c : fixed the angle for bottom one • When two balls collide, flags c is changed to 1, matlab code will jump out of the reducing angle loop

17. Challenge Three -- Tank

18. Throwing Out of the Tank

19. Matlab Code • Compare position of ball with tank wall • Change the velocities when ball hits tank wall • Using flags to switch off the comparison method • If ball is thrown out of the tank, it work as normal, it will hit ground and bound several times

20. ODEs • We are asked to numerically solve (in Matlab) the ordinary differential equation • The exact solution is y(x)=exp(x)-x-1 • We want to solve this problem as quick as possible and also the error is less than 0.1%, therefore Runge Kutta method is used to solve this problem.

21. Matlab Code • Problem • We are asked to find the error at x=10 and keep the error to be less than 0.1%, however, if we want to speed up the code , if 10 is not divisible by h, then we can not obviously show the error at x=10. • Solutions • If the x plus h is greater than 10, which means the next h is out of the range and the difference between x and 10 is equal to 10-x, this is the new step for Runge Kutta. If using the new step to run Runge Kutta, we will calculate y value when x=10.

22. Runge kutta vs. EXACT FUNCTION X increases unit step except last point, when x is close to 10, the increasing step is changed to (10-x)

23. Error for runge kutta method

24. Matlab Code • Step h = 0.3565 • when x=10, • Runge Cutta : y=2.1993e4 • Exact : y=2.2015e4 • error =0.099958% • time take =0.5s