Update 1
My question is this. How do I numerically integrate for the state transition matrix using the differential equations.
Should I be using the Jacobian and inserting data from a reference solution into the jacobian at each time step? My solution of the problem spirals off into infinity, and I've been working on this for a week now. And nobody seems to solve the problem numerically with any examples ever. The only examples I can find are of the time invariant case, which I can solve with standard techniques.
$$ $$
$$ \textbf{Problem Definition } $$
$$ \ddot{r} = – \frac{1}{r^2} \textbf{r} \tag{1} $$
$$ X_o = [x_o \ \ \dot{x_o} \ \ \ y_o \ \ \ \dot{y_o}] = [1 \ \ \ 0 \ \ \ 0 \ \ \ 1] $$
Ignoring the terms
$$ GMm $$
we assume they are unity here.
$$ $$
In component form
$$ \ddot{x} = – \frac{x}{r^3} \tag{2}$$
$$ \ddot{y} = – \frac{y}{r^3} \tag{3}$$
$$ $$
Now the state transition matrix, according to papers I've read, can be defined as below in (4)
$$ \dot\Phi(t,t_o) \ = \ \frac{\partial\dot{x_f}}{\partial{x_f}} \frac{\partial{x_f}}{\partial{x_o}} \tag{4}$$
where,
$$ $$
(According to a few astrodynamics thesis' I've read)
$$ A(t) = \frac{\partial\dot{x_f}}{\partial{x_f}} \tag{5}$$
which is time varying (also the reason I have to integrate the state transition matrix numerically)
$$ $$
And then we have the state transition matrix as,
$$ \Phi(t,t_o) = \frac{\partial{x_f}}{\partial{x_o}} \tag{6} $$
Putting this altogether gives us
$$ $$
$$ \dot\Phi(t,t_o) \ = \ A(t) \ \Phi(t,t_o) \tag{7}$$
My attempt at a solution was to take the partial derivatives of the differential equation with respect to each variable,
$$ x, \ y, \ \dot{x} \ \ \dot{y} $$
making this the jacobian, a matrix of $$ [*]_{4×4} $$
According to what I've read, the jacobian should be evaluated along the 'reference trajectory,' (which is very important, according to books)
What I don't understand is, what value for $$ x_f $$ am I supped to use in evaluating the jacobian?
$$ $$
$$ \textbf{Question } $$
My initial thought was to use the values generated by the true solution, and plug those in for each value of x and y and use that jacobian matrix to integrate for a specific state transition matrix at that one single instant, and update the jacobian at each time step along the way. Is the way to integrate it?
In addition, I'm confused about the notation $$ \dot{x_f}, \ \ and \ \ x_f$$
used in the derivatives of the jacobian and state transition matrix. I really need some clarification on this. How am I supposed to take the derivative of $$\dot{x_f}$$ when I don't really know what $$ x_f $$ is? This goes for $$ x_o $$ as well. In order for me to take the derivative with respect to the initial condition $$ x_o $$ it would have to be present in the equation, but it isn't. Or, am I missing something simple, like, I should just integrate the original differential equation with respect to time (which isn't present) and insert the initial conditions into that boundary value problem and take the derivative.
As you can see, I have several avenues to go from here which I'm confused about.
Best Answer
You need $x$ to be the state (a vector of position and velocity, generally) at some time $t$. Usually when people use the notation $x_f$ it means $x(t_f)$, or the state at time $t_f$.
Then, $\dot{x}_f$ is the vector of instantaneous velocity and acceleration (usually due to gravity at the position given in $x_f$).
So, $A =\frac{\partial \dot{x}_f}{\partial x_f}$ is going to be a Jacobian. For orbital problems it usually looks like $A = \begin{bmatrix}0_{3\times3} & I_{3\times3} \\ G(r) & 0_{3\times3}\end{bmatrix}$, where $r$ is the position component of $x$ at whatever time you're interested in and $G(r)$ is the gravity gradient Jacobian, $\frac{\partial\ddot{r}}{\partial r}$ evaluated at $x$.
You can find a formula for the gravity gradient in a lot of places, but I just happen to have my copy of Battin open (An Introduction to the Mathematics and Methods of Astrodynamics, Revised Edition), and it's given at the bottom of page 454, just under problem 9-11.