- #1
KLoux
- 174
- 1
Hello. I've recently been using a few different numerical integration methods to solve a number of different problems, and as I looked into the integration algorithms, I realized that I don't understand them as well as I thought I did. I created a simple Excel workbook that calculates the motion of a simple pendulum and performed the integration using four different methods: Forward Euler, Second Order Runge Kutta, Third Order Adams Bashforth, and Fourth Order Runge Kutta.
At first I looked only at plots of the pendulum's velocity and position, but all four methods had similar results. I had to look at a very small time scale before I could see the differences. I tried reducing the simulation frequency to see if the differences between the methods would become more obvious. When I dropped the frequency down to 20 Hz, I noticed that that the AB3 simulation was loosing energy. I added a plot to show the energy.
The RK4 method, not surprisingly, appeared to have no energy loss. The RK2 and Euler methods oscillated (Euler had amplitude ~30 times the amplitude of RK2 method) but appeared to maintain a constant average value (no long-term loss). The AB3 method started loosing energy right from the start of the simulation (see attached plot). This brings me to my first questions: Is this a characteristic of the AB methods? Is this something that is specific to this simulation? In another simulation would I be gaining energy?
Initially I thought that the energy loss was primarily dependent on simulation frequency, but I realized that I was simulating a fixed number of time steps rather than a fixed amount of time (a result of using Excel for the simulations). I though perhaps a plot of rate of energy loss would remove this effect. I was surprised to see that the rate of energy loss (of all of the methods) was inversely proportional to the simulation frequency. This makes sense for the methods whose energy oscillates, but I didn't expect it for the RK4 and AB3 methods. I also noticed that the RK4 method was also loosing energy. It was disappearing at a much lower rate (~10e-8 times AB3 - definitely negligible for my purposes) but the rate of change of energy for the method was always negative. So my second questions: Why? I would expect the RK2 and RK4 methods to have similar shapes for these curves, but this is does not appear to be the case. Does it sound like my implementation is flawed?
Most of my simulations run from a C++ application that calls MATLAB to calculate the state derivatives. I previously favored the AB3 method because I can store the state derivatives from previous time steps, and I only need to call MATLAB once every time step to calculate the next time step. This results in a much faster execution and let's me run the simulations at a higher frequency. I am sure that I would have to reduce my simulation frequency to switch to RK4, which would have to make four calls to MATLAB every time step, as I have a requirement to run faster than real-time (and MATLAB is very slow). And now my last questions: Is there a way that I can "add" energy back to the system if I know what my energy loss (due to the integration algorithm) was? Should I look into RK4 or something else that runs faster and will loose less energy? Is energy loss even a valid way to asses the accuracy of a simulation?
If you've read this far through the post, then thank you for your time (and I'm sorry it was so long!). Thanks for your help!
-Kerry
p.s. I was hoping to be able to post the spreadsheet, but the file was too large. If you're interested in seeing it I can remove some plots to make it smaller and try posting it here again, or I can e-mail it to you - just let me know.
At first I looked only at plots of the pendulum's velocity and position, but all four methods had similar results. I had to look at a very small time scale before I could see the differences. I tried reducing the simulation frequency to see if the differences between the methods would become more obvious. When I dropped the frequency down to 20 Hz, I noticed that that the AB3 simulation was loosing energy. I added a plot to show the energy.
The RK4 method, not surprisingly, appeared to have no energy loss. The RK2 and Euler methods oscillated (Euler had amplitude ~30 times the amplitude of RK2 method) but appeared to maintain a constant average value (no long-term loss). The AB3 method started loosing energy right from the start of the simulation (see attached plot). This brings me to my first questions: Is this a characteristic of the AB methods? Is this something that is specific to this simulation? In another simulation would I be gaining energy?
Initially I thought that the energy loss was primarily dependent on simulation frequency, but I realized that I was simulating a fixed number of time steps rather than a fixed amount of time (a result of using Excel for the simulations). I though perhaps a plot of rate of energy loss would remove this effect. I was surprised to see that the rate of energy loss (of all of the methods) was inversely proportional to the simulation frequency. This makes sense for the methods whose energy oscillates, but I didn't expect it for the RK4 and AB3 methods. I also noticed that the RK4 method was also loosing energy. It was disappearing at a much lower rate (~10e-8 times AB3 - definitely negligible for my purposes) but the rate of change of energy for the method was always negative. So my second questions: Why? I would expect the RK2 and RK4 methods to have similar shapes for these curves, but this is does not appear to be the case. Does it sound like my implementation is flawed?
Most of my simulations run from a C++ application that calls MATLAB to calculate the state derivatives. I previously favored the AB3 method because I can store the state derivatives from previous time steps, and I only need to call MATLAB once every time step to calculate the next time step. This results in a much faster execution and let's me run the simulations at a higher frequency. I am sure that I would have to reduce my simulation frequency to switch to RK4, which would have to make four calls to MATLAB every time step, as I have a requirement to run faster than real-time (and MATLAB is very slow). And now my last questions: Is there a way that I can "add" energy back to the system if I know what my energy loss (due to the integration algorithm) was? Should I look into RK4 or something else that runs faster and will loose less energy? Is energy loss even a valid way to asses the accuracy of a simulation?
If you've read this far through the post, then thank you for your time (and I'm sorry it was so long!). Thanks for your help!
-Kerry
p.s. I was hoping to be able to post the spreadsheet, but the file was too large. If you're interested in seeing it I can remove some plots to make it smaller and try posting it here again, or I can e-mail it to you - just let me know.
