Newtonian Motion: Initial-Value and State-Space Formulation

Welcome back. We now have the physical rule needed for the simplest orbital model: in an Earth-centred inertial frame, a satellite at position $\mathbf{r}$ has acceleration

$$\mathbf{a}(\mathbf{r})=-{\mu }_{\mathrm{E}}\frac{\mathbf{r}}{\Vert \mathbf{r}{\Vert }^3}.$$

That rule alone does not yet specify an orbit. To obtain a particular trajectory, we must state where the satellite is and how it is moving at one chosen time, then formulate the problem in the form a numerical solver can evolve.

This lesson turns the second-order Newtonian equation into a first-order state-space initial-value problem. This is the interface behind essentially every ODE solver you will use later: a function takes the present state and returns its instantaneous rate of change.

---

From Newton's law to an initial-value problem

Newton's second law for a satellite of mass $m$ is

$$m\ddot{\mathbf{r}}(t)=\mathbf{F}(t,\mathbf{r},\dot{\mathbf{r}}).$$

For the fixed-Earth point-mass model, the only force is gravity. Dividing by $m$ gives

$$\ddot{\mathbf{r}}(t)=-{\mu }_{\mathrm{E}}\frac{\mathbf{r}(t)}{\Vert \mathbf{r}(t){\Vert }^3}.$$

The double dot means a second derivative with respect to time:

$$\dot{\mathbf{r}}(t)=\mathbf{v}(t),\ddot{\mathbf{r}}(t)=\dot{\mathbf{v}}(t)=\mathbf{a}(t).$$

The equation is second order because its highest time derivative is $\ddot{\mathbf{r}}$. In three dimensions, it represents three coupled scalar equations:

$$\ddot{x}=-{\mu }_{\mathrm{E}}\frac{x}{(x^2+y^2+z^2{)}^{3/2}},$$ $$\ddot{y}=-{\mu }_{\mathrm{E}}\frac{y}{(x^2+y^2+z^2{)}^{3/2}},$$ $$\ddot{z}=-{\mu }_{\mathrm{E}}\frac{z}{(x^2+y^2+z^2{)}^{3/2}}.$$

These are nonlinear equations: the right-hand side is not a fixed matrix times the coordinate vector, because the denominator depends nonlinearly on the present position.

A differential equation describes which curves are physically allowed. An initial-value problem selects one particular allowed curve by declaring the state at an initial time $t_0$:

$$\mathbf{r}(t_0)={\mathbf{r}}_0,\dot{\mathbf{r}}(t_0)={\mathbf{v}}_0.$$

Thus the complete restricted two-body initial-value problem is

$$\boxed{\begin{array}{rl}\ddot{\mathbf{r}}(t)& =-{\mu }_{\mathrm{E}}\frac{\mathbf{r}(t)}{\Vert \mathbf{r}(t){\Vert }^3},\\ \mathbf{r}(t_0)& ={\mathbf{r}}_0,\\ \dot{\mathbf{r}}(t_0)& ={\mathbf{v}}_0.\end{array}}$$

The initial position supplies three scalar values, and the initial velocity supplies another three. That is exactly the amount of initial information required for three second-order equations.

It is worth distinguishing this from a boundary-value problem. In an initial-value problem, all known conditions occur at one time, such as launch time. In a boundary-value problem, one might instead prescribe a position at launch and another position at a later arrival time. Orbital propagation is naturally an initial-value problem: given a state estimate now, predict its future evolution.

---

What counts as the state?

A system's state is the smallest collection of quantities that, together with the governing model, determines its instantaneous future evolution.

For the satellite, position alone is not enough. Two spacecraft can occupy exactly the same point in space while travelling with different velocities, and their subsequent paths will differ. Acceleration is not an independent part of the state in this model, because it is already determined by position through the gravitational law.

Define the six-dimensional state vector

$$\mathbf{y}(t)=\left[\begin{array}{c}\mathbf{r}(t)\\ \mathbf{v}(t)\end{array}\right]=\left[\begin{array}{c}x(t)\\ y(t)\\ z(t)\\ v_x(t)\\ v_y(t)\\ v_z(t)\end{array}\right].$$

The set of possible states is called state space, or, in mechanics, phase space. Here it is effectively six-dimensional position-velocity space, excluding the singular state $\mathbf{r}=0$, where the point-mass gravitational acceleration is undefined.

At time $t_0$, the initial state is simply

$$\boxed{\mathbf{y}(t_0)={\mathbf{y}}_0=\left[\begin{array}{c}{\mathbf{r}}_0\\ {\mathbf{v}}_0\end{array}\right].}$$

For example, using kilometres and seconds, a possible initial state might be assembled as

$${\mathbf{r}}_0=\left[\begin{array}{c}7000\\ 0\\ 0\end{array}\right]\mathrm{km},{\mathbf{v}}_0=\left[\begin{array}{c}0\\ 7.5\\ 0\end{array}\right]\mathrm{km}{\mathrm{s}}^{-1}.$$

This specifies a satellite initially $7000$ km from Earth's centre on the positive $x$-axis, travelling initially in the positive $y$-direction. At this stage, regard it as an input state rather than trying to classify its resulting orbit. Later, you will derive how position and velocity determine such features as eccentricity and semi-major axis.

---

Converting the second-order equation to first-order state space

The conversion is conceptually simple but foundational. Introduce velocity as a separate state variable:

$$\mathbf{v}(t)=\dot{\mathbf{r}}(t).$$

This gives one first-order equation:

$$\dot{\mathbf{r}}=\mathbf{v}.$$

The gravitational law gives the other:

$$\dot{\mathbf{v}}=-{\mu }_{\mathrm{E}}\frac{\mathbf{r}}{\Vert \mathbf{r}{\Vert }^3}.$$

Stacking these equations yields

$$\boxed{\frac{d}{dt}\left[\begin{array}{c}\mathbf{r}\\ \mathbf{v}\end{array}\right]=\left[\begin{array}{c}\mathbf{v}\\ -{\mu }_{\mathrm{E}}\mathbf{r}/\Vert \mathbf{r}{\Vert }^3\end{array}\right].}$$

Equivalently, in the standard numerical-analysis notation,

$$\boxed{\dot{\mathbf{y}}=f(t,\mathbf{y}),\mathbf{y}(t_0)={\mathbf{y}}_0,}$$

where, for this autonomous gravitational model,

$$f(t,\mathbf{y})=f(\mathbf{y})=\left[\begin{array}{c}\mathbf{v}\\ -{\mu }_{\mathrm{E}}\mathbf{r}/\Vert \mathbf{r}{\Vert }^3\end{array}\right].$$

The function is called autonomous because its formula does not depend explicitly on time. It depends only on the current state. A solver conventionally still passes $t$ as an argument because many physically important models do depend explicitly on time: thrust schedules, time-varying forcing, or a prescribed rotating frame are examples.

In scalar components, the same first-order system is

$$\frac{d}{dt}\left[\begin{array}{c}x\\ y\\ z\\ v_x\\ v_y\\ v_z\end{array}\right]=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -{\mu }_{\mathrm{E}}x/(x^2+y^2+z^2{)}^{3/2}\\ -{\mu }_{\mathrm{E}}y/(x^2+y^2+z^2{)}^{3/2}\\ -{\mu }_{\mathrm{E}}z/(x^2+y^2+z^2{)}^{3/2}\end{array}\right].$$

Notice the structure:

Entries of $\mathbf{y}$	Their time derivatives
Position $\mathbf{r}$	Velocity $\mathbf{v}$
Velocity $\mathbf{v}$	Acceleration $\mathbf{a}(\mathbf{r})$

That is the whole reduction. A second-order problem has not been approximated or simplified physically; it has been rewritten exactly in a form that exposes all quantities needed to continue the motion from the current state.

More generally, any second-order vector equation

$$\ddot{\mathbf{q}}=\mathbf{F}(t,\mathbf{q},\dot{\mathbf{q}})$$

can be written using $\mathbf{u}=\dot{\mathbf{q}}$ as

$$\frac{d}{dt}\left[\begin{array}{c}\mathbf{q}\\ \mathbf{u}\end{array}\right]=\left[\begin{array}{c}\mathbf{u}\\ \mathbf{F}(t,\mathbf{q},\mathbf{u})\end{array}\right].$$

This general pattern is part of the general-purpose applied-mathematics toolkit: it applies to mechanical systems, coupled oscillators, electrical circuits expressed through suitable state variables, and many nonlinear dynamical systems.

---

The reading's twelve-component formulation is not a competing notation for the current model. It is the model you will reach in the next module when Earth is no longer held fixed. The state-space architecture already scales cleanly:

• fixed central body and one satellite: $6$ state components;
• two freely moving bodies: $12$ state components;
• $N$ freely moving point masses: $6N$ state components.

Only the force calculation changes substantially.

---

A compact Python representation

Your right-hand-side function should implement the mathematical object $f(t,\mathbf{y})$. It does not advance time itself. It answers a narrower question:

Given the current state, what is its instantaneous derivative?

import numpy as np

def earth_orbit_rhs(t, state, mu):
    """
    Return d(state)/dt for a test satellite orbiting a fixed point mass.

    state = [x, y, z, vx, vy, vz]
    position units: km
    velocity units: km / s
    mu units: km**3 / s**2
    """
    position = state[:3]
    velocity = state[3:]

    r2 = np.dot(position, position)

    if r2 == 0.0:
        raise ValueError("Point-mass gravitational acceleration is undefined at r = 0.")

    acceleration = -mu * position / (r2 * np.sqrt(r2))

    return np.concatenate((velocity, acceleration))

The returned array is

$$\dot{\mathbf{y}}=\left[\begin{array}{c}\mathbf{v}\\ \mathbf{a}\end{array}\right],$$

not the next position and velocity. This distinction is crucial. An ODE integrator will later call this function repeatedly, combining derivative evaluations according to a particular numerical method.

You can construct an initial state in the same ordering:

r0 = np.array([7000.0, 0.0, 0.0])  # km
v0 = np.array([0.0, 7.5, 0.0])     # km / s

state0 = np.concatenate((r0, v0))

A useful unit audit of the return value is:

$$[\dot{\mathbf{y}}]=\left[\frac{\mathrm{km}}{\mathrm{s}},\frac{\mathrm{km}}{\mathrm{s}},\frac{\mathrm{km}}{\mathrm{s}},\frac{\mathrm{km}}{{\mathrm{s}}^2},\frac{\mathrm{km}}{{\mathrm{s}}^2},\frac{\mathrm{km}}{{\mathrm{s}}^2}\right].$$

It is normal for a state vector to contain quantities with different units. Its derivative correspondingly contains velocity units for the position entries and acceleration units for the velocity entries.

Three implementation errors are especially common:

1. Reversing the blocks. Returning $[\mathbf{a},\mathbf{v}]$ is wrong for the declared ordering $[\mathbf{r},\mathbf{v}]$. The first three derivatives must be velocities.

2. Using position alone as the initial state. A position gives the gravitational acceleration, but it does not determine the orbit. The velocity vector is equally necessary.

3. Treating the state equation as linear. A linear system can be written as $\dot{\mathbf{y}}=A\mathbf{y}$ for a constant matrix $A$. Gravity cannot: its acceleration coefficient changes with $\Vert \mathbf{r}\Vert$.

---

---

Why this form is the computational boundary

The state-space equation establishes a clean separation between three tasks:

1. The physical model determines $f$. Here, it is the inverse-square gravitational acceleration.

2. The initialisation determines ${\mathbf{y}}_0$. Here, it is an initial position and velocity in a clearly specified inertial frame and unit system.

3. The numerical method estimates the trajectory $\mathbf{y}(t)$ from the initial-value problem. We have deliberately not chosen that method yet.

This separation is valuable for scrutiny. If an orbit looks wrong, you can ask a more precise question: was the force law coded incorrectly, were the initial conditions inconsistent, or did the integrator introduce excessive numerical error? Those are distinct failure modes, and a well-structured simulator makes them distinguishable.

There is also a mathematical domain condition. The right-hand side is smooth for every state with

$$\Vert \mathbf{r}\Vert >0,$$

but it becomes singular at the origin. Given a valid initial state away from that singularity, the deterministic model specifies the local future rate of change uniquely. A physical Earth-orbit simulation should additionally stop or change its model before a trajectory reaches Earth's surface; the point-mass equation alone is not a collision or atmospheric-entry model.

---

Key takeaways

Newtonian orbital motion in the fixed-Earth model begins as a second-order vector equation:

$$\ddot{\mathbf{r}}=-{\mu }_{\mathrm{E}}\frac{\mathbf{r}}{\Vert \mathbf{r}{\Vert }^3}.$$

A specific orbit requires six initial scalar values:

$$\mathbf{r}(t_0)={\mathbf{r}}_0,\mathbf{v}(t_0)={\mathbf{v}}_0.$$

Defining the state vector as

$$\mathbf{y}=\left[\begin{array}{c}\mathbf{r}\\ \mathbf{v}\end{array}\right]$$

converts the equation exactly into the first-order state-space form

$$\dot{\mathbf{y}}=\left[\begin{array}{c}\mathbf{v}\\ -{\mu }_{\mathrm{E}}\mathbf{r}/\Vert \mathbf{r}{\Vert }^3\end{array}\right],\mathbf{y}(t_0)={\mathbf{y}}_0.$$

This right-hand-side function is the central software boundary for your simulator: it maps a current state to its instantaneous derivative, while a later integrator will decide how to use those derivatives to move through time.

Next, we will nondimensionalise the gravitational model. That will reveal its natural orbital timescale and give you a cleaner way to reason about parameter magnitudes, timestep choices, and simulations of systems far beyond Earth.