14.5: Conservation of Energy

Last updated
Save as PDF

Page ID: 84758

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

A useful way to check the accuracy of an ODE solver is to see whether it conserves energy. For planetary motion, it turns out that ode45 does not.

The kinetic energy of a moving body is

\[KE = m v^2 / 2\notag \]

The potential energy of a sun with mass \(m_1\) and a planet with mass \(m_2\) and a distance \(r\) between them is

\[PE = -G \frac{m_1 m_2}{r}\notag\]

Exercise 14.5

Write a function called energy_func that takes the output of your Earth simulation and computes the total energy (kinetic and potential) of the system for each estimated position and velocity.

Plot the result as a function of time and check whether it increases or decreases over the course of the simulation.

You can reduce the rate of energy loss by decreasing ode45’s tolerance option using odeset (see “ODE Events” on page 11.3):

options = odeset('RelTol', 1e-5);
[T, M] = ode45(@rate_func, tspan, W, options);

The name of the option is RelTol for “relative tolerance.” The default value is 1e-3, or 0.001. Smaller values make ode45 less “tolerant,” so uses smaller step sizes to make the errors smaller.

Run ode45 with a range of values for RelTol and confirm that as the tolerance gets smaller, the rate of energy loss decreases.

Along with ode45, MATLAB provides several other ODE solvers (see https://greenteapress.com/matlab/solver). Run your simulation with one of the other ODE solvers MATLAB provides and see if any of them conserve energy. You might find that ode23 works surprisingly well (although technically it does not conserve energy either).