Problem-solving tips¶

Many problems we face will require common solution approaches involving solving equations (i.e., finding the root of an equation) and integrating systems of ordinary differential equations (ODEs). This module shows some tips and strategies, and examples, for doing this.

Solving equations numerically¶

In many problems we need to solve an equation. For example, let’s think about solving \(f(x) = g\). This means we want to find the value \(x\) that satisfies this equation. One way to to think about this is that we want find the root of \(f(x) - g = 0\).

When the equation is complicated, nonlinear, and/or transcendental, we probably cannot find the root analytically, and so we need to solve for it numerically.

In Matlab, we can use the fzero() function, while in Python we can use the root_scalar() function provided in SciPy’s optimize module.

For example, let’s see how to solve the equation for eccentric anomaly \(E\), part of the Kepler equations for an elliptic orbit:

\[ M_e = 2 \pi \frac{t}{T} = E - e \sin E \;, \]

where \(t\) is the time of flight, \(T\) is the orbital period, \(e\) is the orbital eccentricity, and \(M_e\) is the mean anomaly.

To solve any nonlinear equation like this using a numerical method, we need to shift the equation to a zero-valued form, like \(f(x) = 0\), where \(x\) is the root of the equation and what we are trying to find. In this case, that means the equation looks like

\[ 2 \pi \frac{t}{T} - \left( E - e \sin E \right) = 0 \;. \]

We also need to identify a reasonable first guess for the numerical method to start its iterations. In this case, that is the Prussing estimate: \(E_0 = M_e + \frac{e}{2}\). But, in general this first guess should be a value that is reasonable based on the problem at hand.

Now, let’s look at how to solve this:

matlab

We need to define the known values to calculate the initial guess, and also make sure the necessary values are in the function.

e = 0.1;
t = 10;
T = 60;
M_e = 2 * pi * t / T;
E0 = M_e + e/2;

E = fzero(@f, E0);

function zero = f(E)
    e = 0.1; t = 10; T = 60;
    zero = (2*pi*t/T) - (E - e * sin(E));
end

Note that this code snippet assumes the function is defined inside your script file, at the bottom. You could also define it in a separate .m file, but then you would need to call it using E = fzero('name.m', E0);.

This approach does rely on us defining variables twice, which is a pain if we need to change any of them. We can be a bit clever and instead define a new parameterized function:

e = 0.1;
t = 10;
T = 60;
M_e = 2 * pi * t / T;
E0 = M_e + e/2;

f = @(E) eccentric_anomaly(E, e, t, T);
E = fzero(f, E0);

function zero = eccentric_anomaly(E, e, t, T)
    zero = (2*pi*t/T) - (E - e * sin(E));
end

python

We can define the constants, and then use those both inside the function and in the main scope.

import numpy as np
from scipy.optimize import root_scalar

# given values
e = 0.1
t = 10
T = 60

def eccentric_anomaly(E):
   return (2*np.pi*t/T) - (E - e * np.sin(E))

M_e = 2 * pi * t / T
E0 = M_e + e/2
sol = root_scalar(eccentric_anomaly, x0=E0)

We can also be more careful and include the necessary variables at parameters to the function:

import numpy as np
from scipy.optimize import root_scalar

def eccentric_anomaly(E, e, t, T):
   return (2*np.pi*t/T) - (E - e * np.sin(E))

# given values
e = 0.1
t = 10
T = 60

M_e = 2 * pi * t / T
E0 = M_e + e/2
sol = root_scalar(eccentric_anomaly, x0=E0, args=(e, t, T))

Space Systems Engineering notes

Problem-solving tips¶

Solving equations numerically¶