# Numerical integration with variation of parameters

I'm having a lot of trouble with a numerical integration problem for an astrodynamics course. We're starting with the function for a perturbed oscillator

$$\ddot{x} + 3 x + x^3 = 0$$

whose homogenous counterpart is

$$\ddot{x} + 3 x = 0$$

and we're told to treat $$k_3x^3$$ as the perturbation. The goal is to set up a variation of parameters problem and solve for the parameters themselves (to determine what value they would have to be at a certain time in order to hold the unperturbed solution's form fixed). I believe the unperturbed/homogenous solution is:

$$x(t) = c_1 \cos(\sqrt{3}t)+c_2\sin(\sqrt{3}t)$$

The parameter $$c_1$$ and $$c_2$$ are set to the time varying parameter such that:

$$e_1(t) = c_1$$ $$e_2(t) = c_2$$

and we're given initial conditions

$$e_1(t) = 1$$ $$e_2(t) = 0$$

So I want to find expressions for these time varying parameters such that

$$x(t) = e_1(t) \cos(\sqrt{3}t)+e_2(t)\sin(\sqrt{3}t)$$

holds. I took the derivative of this and set up the osculating conditions (the constraints on the varying parameters to match the perturbed and unperturbed expressions) to get

$$\dot{x}(t) = e_1[-\sqrt{3}\sin(\sqrt{3}t)]+e_2[\sqrt{3}\cos(\sqrt{3}t)]$$

I've tried to figure this out using the Wronskian, but I'm a bit confused since on the RHS we have $$0$$, and since our perturbation has an $$x$$ in it (and will vary with time)? Following the method from our class (or just computing the entire first and second derivatives), I get the constraints

$$\dot{e}_1\cos(\sqrt{3}t)+\dot{e}_2\sin(\sqrt{3}t) = 0$$ $$\dot{e}_1[-\sqrt{3}\sin(\sqrt{3}t)]+\dot{e}_2[\sqrt{3}\cos(\sqrt{3}t)] = -x^3$$

which we were taught to solve by setting up a matrix equation, but that I've also tried by solving for one and substituting. $$\begin{bmatrix} \cos(\sqrt{3}t) & \sin(\sqrt{3}t) \\ -\sqrt{3}\sin(\sqrt{3}t) & \sqrt{3}\cos(\sqrt{3}t) \end{bmatrix} \begin{bmatrix} \dot{e}_1 \\ \dot{e}_2 \end{bmatrix} = \begin{bmatrix} 0 \\ -x^3 \end{bmatrix}$$

Inverting:$$\begin{bmatrix} \dot{e}_1 \\ \dot{e}_2 \end{bmatrix} = \begin{bmatrix} (\cos(\sqrt{3} t) & \frac{-\sin(\sqrt{3} t)}{\sqrt{3}} \\ \sin(\sqrt{3} t) & \frac{\cos(\sqrt{3} t)}{\sqrt{3}} \end{bmatrix} \begin{bmatrix} 0 \\ -x^3 \end{bmatrix}$$

Which gives me the differential equations

$$\dot{e}_1 = \frac{\sin(\sqrt{3} t) x^3}{\sqrt{3}}$$ $$\dot{e}_2 = -\frac{\cos(\sqrt{3} t) x^3}{\sqrt{3}}$$

And this is where I get confused (assuming I've done the previous work correctly). The $$x$$ is that expression from earlier that depends on time. I need to numerically integrate these and get there value after 10 seconds, but $$x(t)$$ already has $$e_1$$ and $$e_2$$ in it

$$x(t) = e_1(t) \cos(\sqrt{3}t)+e_2(t)\sin(\sqrt{3}t)$$

which is what we're solving for. Seeing this, I imagine that we must instead be using the unperturbed expression for $$x(t)$$ that has the initial conditions plugged in, but I'm still getting the wrong answer. I'm conceptually a little on clear on how to integrate these. I've been using python's solve_ivp, and as far as I know these can be integrated separately with their functional form and the provided initial conditions. I've done tons of variation of parameters practice problems to try to help with this which I ran into no trouble with. Any advice?

Here's my code: interestingly, I've gotten the first parameter right, but am somehow wrong on the second.

Python code

Here is an image of the question itself: original question

I’ve also plotted the solution for the homogenous case against the perturbed one: plot

• Do question and answer in math.stackexchange.com/questions/3961162/…, math.stackexchange.com/questions/2013417/… help you along? The main change in thinking is that you make a series expansion $x=x_0+k_3x_1+...$ with $k_3$ as very small parameter Commented Apr 11 at 7:40
• @LutzLehmann I appreciate the response and the links!! Btw, I uploaded an image of the code I'm using to integrate this: I get the first parameter correct and the second incorrect, sadly (which is weird, as I'm not sure how I get one parameter without the other). I might be misunderstanding and I'm rusty on perturbative expansions, but would not expanding as a series in powers of $k_3$ just be an alternative to the variation of parameters method I'm using here? I guess I would expand in $k_3$ rather than t, and match coefficients? Thanks again, btw. Commented Apr 11 at 16:04
• No, that is separate. The separation approach gives you a series of differential equation with a left side that is always identical and a right side that only contains solutions to previous equations of the sequence. While more general, using variation is tedious. In examples that can be computed manually like this one, most often the guessing method, "undetermined parameters" is applicable Commented Apr 11 at 16:09
• @LutzLehmann ah yeah, it’s not my first choice. In this case though what the prof wants are the values of the parameters that we’re varying to keep the form of the perturbed solutions the same as that of unperturbed ones—rather than the solution to the perturbed problem itself that you’d get in the end. Is that something you can get from series? Commented Apr 11 at 16:21
• So the aim is really that you want to get an alternative first-order system that is connected to a co-rotating coordinate system. A transformation to polar coordinates would go one step beyond that. But then you do not need a perturbation idea. But as perturbation is mentioned in the task, it hints that some kind of perturbation is expected in the solution. Could you edit your question to sharpen andput in full detail what the task demanded and what result is expected? Your theory looks good as written. Commented Apr 11 at 17:28

## 1 Answer

Using the Lindsted-Poincaré method, from $$\ddot u+\omega_0^2u + \epsilon u^3=0$$ making $$\omega t = \tau$$ we have

$$\omega^2u''+\omega_0u+\epsilon u^3=0$$

considering now $$\omega = \omega_0 + \epsilon \omega_1 + \epsilon^2 \omega_2+\cdots$$ as well as $$u = u_0 + \epsilon u_1 + \epsilon^2 u_2+\cdots$$ after substitution we have

$$\left(\sum_{k=0}^n\epsilon^k\omega_k\right)^2\left(\sum_{j=0}^m \epsilon^k u_k\right)''+\omega_0\left(\sum_{j=0}^n \epsilon^k u_k\right)+\epsilon\left(\sum_{j=0}^n \epsilon^k u_k\right)^3=0$$

now considering $$n=1$$ and developing in powers of $$\epsilon$$ we have

$$\cases{ \ddot u_0 + u_0 = 0\\ 2 \omega _0 \omega _1 \ddot u_0+\omega _0^2 \ddot u_1+\omega _0^2 u_1+u_0^3 = 0 }$$

from the first ode we have $$u_0 = a_0\cos(\tau+b_0)$$ and substituting into the second ode

$$\omega _0^2 \ddot u_1+\omega _0^2 u_1+a_0^3 \cos ^3(b_0+\tau)-2 a_0 \omega _0 \omega _1 \cos (b_0+\tau)=0$$

and solving this ode we have

$$u_1 = a_1\cos(\tau+b_1) -\frac{3 a_0^3 \tau \sin (b_0+\tau)}{8 \omega _0^2}+\frac{a_0 \omega _1 \tau \sin (b_0+\tau)}{\omega _0}-\frac{3 a_0^3 \cos (b_0+\tau)}{16 \omega _0^2}+\frac{a_0^3\cos (3 (b_0+\tau))}{32 \omega _0^2}+\frac{a_0 \omega _1 \cos (b_0+\tau)}{2 \omega _0}$$

now we can eliminate the secular terms by making

$$\omega_1 = \frac{3a_0^2}{8\omega_0}$$

resulting in

$$u_1 = a_1\cos(\tau+b_1) +\frac{a_0^3}{32\omega_0^2}\cos(3(\tau+b_0))$$

so we got with $$n=1$$

$$u = u_0 + \epsilon u_1$$

Note that $$\tau = \omega t = (\omega_0+\epsilon \omega_1)t$$. Any way, the solutions are periodic because from $$\ddot u+\omega_0^2u + \epsilon u^3=0$$ after multiplication by $$\dot u$$ and integration, we have

$$\frac 12 \dot u^2+\frac 12\omega_0^2u^2+\frac{\epsilon}{4}u^4 = c_0$$

which describes periodic orbits on the phase plane.