1 / 38

Numerical Methods for Ordinary Differential Equations and Gravitational Simulations

This lecture reviews numerical schemes for solving ordinary differential equations (ODEs), with a focus on the Forward Euler method and its application in Newtonian motion and gravitation. Key concepts include stability, convergence, and the use of MATLAB for simulating gravitational interactions among particles. Attendees will learn about different numerical methods like Runge-Kutta and Adams-Bashforth schemes, error estimation, and practical applications involving planetary motion. The session will also involve a team exercise to run a planetary simulation script.

rey
Télécharger la présentation

Numerical Methods for Ordinary Differential Equations and Gravitational Simulations

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. MA/CS 375 Fall 2002 Lecture 13 MA/CS 375 Fall 2002

  2. Office Hours • My office hours • Rm 435, Humanities, Tuesdays from 1:30pm to 3:00pm • Rm 435, Humanities, Thursdays from 1:30pm to 3:00pm • Tom Hunt is the TA for this class. His lab hours are now as follow • SCSI 1004, Tuesdays from 3:30 until 4:45 • SCSI 1004, Wednesdays from 12:00 until 12:50 • Hum 346 on Wednesdays from 2:30-3:30 MA/CS 375 Fall 2002

  3. Review Ordinary Differential Equation • Example: • we should know from intro calculus that: • then obviously: MA/CS 375 Fall 2002

  4. Review Forward Euler Numerical Scheme • Numerical scheme: • Discrete scheme: where: MA/CS 375 Fall 2002

  5. Review Direct Proof of Convergence • Fix T and let delta t drop to zero • In this case n needs to increase to infinity MA/CS 375 Fall 2002

  6. Review Stable Approximations • 0<dt<1 dt dt=0.125 dt=0.5 dt=0.25 MA/CS 375 Fall 2002

  7. Review Application: Newtonian Motion MA/CS 375 Fall 2002

  8. m2 m1 Each particle has a scalar mass quantity Review Two Gravitating Particle Masses MA/CS 375 Fall 2002

  9. Review Particle Positions Each particle has a vector position x2 x1 (0,0) MA/CS 375 Fall 2002

  10. Review Particle Velocities Each particle has a vector velocity v1 v2 MA/CS 375 Fall 2002

  11. Review Particle Accelerations Each particle has a vector acceleration a1 a2 MA/CS 375 Fall 2002

  12. Review N-Body Newtonian Gravitation • For particle n out of N The force on each particle is a sum of the gravitational force between each other particle MA/CS 375 Fall 2002

  13. Review N-Body Newtonian Gravitation Simulation • Goal: to find out where all the objects are after a time T • We need to specify the initial velocity and positions of the objects. • Next we need a numerical scheme to advance the equations in time. • Can use forward Euler…. as a first approach. MA/CS 375 Fall 2002

  14. Review Numerical Scheme For m=1 to FinalTime/dt For n=1 to number of objects End For n=1 to number of objects End End MA/CS 375 Fall 2002

  15. Review planets1.m Matlab script • I have written a planets1.m script. • The quantities in the file are in units of • kg (kilograms -- mass) • m (meters – length) • s (seconds – time) • It evolves the planet positions in time according to Newton’s law of gravitation. • It uses Euler-Forward to discretize the motion. • All planets are lined up at y=0 at t=0 • All planets are set to travel in the y-direction at t=0 MA/CS 375 Fall 2002

  16. Mercury has nearly completed its orbit. Data shows 88 days. Run for 3 more days and the simulation agrees!!!. Earth Venus Review Sun Mercury MA/CS 375 Fall 2002

  17. Today Team Exercise • Get the planets1.m file from the web site • This scripts includes: • the mass of all planets and the sun • their mean distance from the sun • the mean velocity of the planets. • Run the script, see how the planets run! • Add a comet to the system (increase Nsphere etc.) • Start the comet out near Jupiter with an initial velocity heading in system. • Add a moon near the earth. • Extra credit if you can make the comet loop the sun and hit a planet  MA/CS 375 Fall 2002

  18. This Lecture • More accurate schemes • More complicated ODEs • Variable time step and embedded methods used to make sure errors are within a tolerance. MA/CS 375 Fall 2002

  19. Adams-Bashforth Schemes • In the forward Euler scheme we only used the value of the right hand side at the previous time step. • i.e. we only used a linear approximation to the time derivative MA/CS 375 Fall 2002

  20. AB Schemes • Idea: we set • where If interpolates fn,fn-1,fn-2,..,fn+1-Nstages • i.e.: MA/CS 375 Fall 2002

  21. f0 f1 f2 f3 We interpolate the function through thefirst 4 points. Then we integrate under the curve between t=3 and t=4. MA/CS 375 Fall 2002

  22. AB Schemes Essentially we use interpolation and a Newton-Cotes quadrature formula to formulate:

  23. Runge-Kutta Schemes • See van Loan for derivation of RK2 and RK4. • I prefer the following (simple) scheme due to Jameson, Schmidt and Turkel (1981): MA/CS 375 Fall 2002

  24. Runge-Kutta Schemes • Beware, it only works when f is a function of y and not t • here s is the order of the scheme. MA/CS 375 Fall 2002

  25. Error Estimate • Matlab has a number of time integrators built in. • ode23 • ode45 • and others.. MA/CS 375 Fall 2002

  26. ode23 • For n=1:#timesteps • ode23 uses two estimates for yn+1 . • A 2nd order RK scheme and a 3rd order RK scheme are used to build two guesses for yn+1. • If the difference between these two estimates are within a tolerance ode23 progresses on to calculating yn+2 • If the difference is greater than the specified tolerance, ode23 reduces the dt and tries again. It repeats until the difference is lower than the tolerance. End MA/CS 375 Fall 2002

  27. Planets Example Using ode23 • Idea: replace our home grown Euler Forward scheme with Matlab’s ode23 scheme in the planets1.m script. MA/CS 375 Fall 2002

  28. Parameters forthe planets asbefore MA/CS 375 Fall 2002

  29. Initial velocities Gather X,Y,VX,VY into one vector MA/CS 375 Fall 2002

  30. Call to Matlab ode23 function [T,CoordVel] = ode23(‘forcing’, TimeLimits, CoordVel(:)); function whichcalculates f(y) Vector holds X,Y,VX,VY [0 FinalTime] Works in Matlab v6.1.0 … at least MA/CS 375 Fall 2002

  31. Call to ode23 Find out the time at each time step Extract final coordinates and velocities Plot planet orbits MA/CS 375 Fall 2002

  32. The routine whichcalculates the forcing function. MA/CS 375 Fall 2002

  33. Team Exercise • Grab planets2.m and forcing.m from the http://useme.org/MA375.html • Run the script • Use the Tsteps vector to find out the time step for each integration stage. Plot a graph showing the time step (dt) at each time step. • Use help to find out how to change the tolerance used by ode23 (hint you will need to use odeset) • Rerun the simulation with a tolerance of 0.1 MA/CS 375 Fall 2002

  34. Application: One-Dimensional Electrostatic Motion MA/CS 375 Fall 2002

  35. Charge Repulsion • Now we will consider the case of charged particles with the same sign charge • Instead of attracting each other, the charges repel each other. MA/CS 375 Fall 2002

  36. Particle Accelerations Each particle has a vector accelerationdirectly away from the other particle a2 a1 MA/CS 375 Fall 2002

  37. Team Project Q1) Modify the planets2.m and forcing.m to simulate the following: • There are N electrically charged particles confined to move in the x-direction only • Distribute the charges initially at equispaced points in [-1,1] • The equations of motion of the charges are: MA/CS 375 Fall 2002

  38. Team ProjectReport Required on 9/25/02 Q1cont) Include the following in one project write up per team: • A scatter plot, with x as the horizontal axis and t as the vertical axis, showing the paths of all the charged particles using ode23 • A graph showing the size of each time step used by ode23 • Replace ode23 with ode45 and rerun • Plot the ode23 and ode45 (t,x) paths on the same graph. • Plot the ode23 and ode45 time step sizes on the same graph • Names of team members MA/CS 375 Fall 2002

More Related