# Homogeneous linear differential equation system with simple periodical coefficient matrix

Hello, I encountered the following system of linear first-order differential equations:

$y'(z)=A(z) y(z)$

where

$y(z): R \rightarrow R^2$ and $A(z)=\begin{pmatrix} 0 & B Cos(\alpha z + \Phi_b) \cr A Cos(\alpha z + \Phi_a) & 0 \end{pmatrix}$

All quantities are real and constant if not explicitly variable. (Probably - if anyone wants to talk about approximations, which would be ok to some extent -) $A$,$B$ are $<1$, while $\alpha>10^3$ and $z<10^{-2}$.

I found Floquet's theorem, which only gives my the structure of the solution but isn't of much help for me. In approximation, I thought about using the Magnus expansion.

I added math.SP, as I got the suggestion to transform the equation to a Schrodinger equation with $Y(z)=A_{1,2}(z)^{1/2} y(z)$ but that did not produce any new ideas for a solution.

I am at the end my (physicist's) knowledge, hope anyone has an idea.

-
Does $\Phi_a$ need to be different form $\Phi_b$? If not, then for rescaling DSolve[{y1'[t] == b Cos[t] y2[t], y2'[t] == a Cos[t] y1[t]}, {y1[t], y2[t]}, t] Mathematica returns {{y1[t] -> C[1] Cosh[Sqrt[a] Sqrt[b] Sin[t]] + I C[2] Sinh[Sqrt[a] Sqrt[b] Sin[t]], y2[t] -> (I Sqrt[a] C[2] Cosh[Sqrt[a] Sqrt[b] Sin[t]])/Sqrt[b] + ( Sqrt[a] C[1] Sinh[Sqrt[a] Sqrt[b] Sin[t]])/Sqrt[b]}} –  Piotr Migdal May 27 '11 at 12:16
If $\Phi_a=\Phi_b$, we can transform to a constant coefficient system by a change of the independent variable. Diagonalizing $A(z)$ does NOT decouple the system! –  Michael Renardy May 27 '11 at 12:28
The monodromy matrix in this case is self-adjoint (up to conjugation), so you'll get one exponentially growing solution and one exponentially decaying one with slight periodic oscillations in each. Does this make sense to you? If it does, I'll post the details (how to approximate exponents, oscillation parameters, eigenvectors, etc.). If it doesn't, you'd better reconsider your model. –  fedja May 27 '11 at 18:07
@fedja: Yes, that does make sense. As those are amplitudes, I can make a guess for their squares, the intensities. For the boundary values $y(0)=(1,0)$ (which are the ones I need, I should have added that) and $a=0$, the problem is simple and gives just an oszillation, which is ok. For $a\neq0$, I expect and exponential-like increase in the second element of $y(t)$ and then a small oscillation around a fixed value. Of cause, I did some numerical tests and it seems to work. So yes, your solutions sounds promising. –  elcron May 28 '11 at 6:37
@ Piotr and Michael Yes, I fear those have to be different. It is a simplification from a slightly more complicated matrix where the elements are $A Cos(a t + \Phi_A) + B Sin(a t + \Phi_B)$ and $A Cos(a t + \Phi_A) - B Sin(a t + \Phi_B)$. Here, we can talk about $\Phi_A=\Phi_B$ but definitively $A\neq B$. (As I don't have the worksheet until monday, I am reconstructing this from my head, it might be slightly different, so don't put to much efford in that.) –  elcron May 28 '11 at 6:37
add comment

## 3 Answers

OK. First of all, change $(y_1,y_2)$ to $(y_1,\sqrt{B/A}y_2)$ and the time $z$ to $t=\alpha z+\frac{\Phi_A+\Phi_B}{2}$. Then we'll get the system with the matrix $$A(t)=\begin{pmatrix} 0 & C\cos (t + \Phi) \cr C \cos(t- \Phi) & 0 \end{pmatrix}$$ where $C=\alpha^{-1}\sqrt{AB}$ and $\Phi=(\Phi_B-\Phi_A)/2$. We want our approximation to be decent on $[\alpha z_{\min} ,\alpha z_{\max}]\subset [-\alpha,\alpha]$. Noting that $\cos(t-\Phi)=\cos(\Phi-t)$, we see that $A(t)=A^*(2\pi-t)$, which immediately tells us that the monodromy matrix $M$ from $0$ to $2\pi$ is self-adjoint. Also, denoting $\psi(t)=C\cos(\Phi+t)$=ce^{it}+\bar c{e^{-it}}$with$c=\frac 12Ce^{i\Phi}, we see that the fundamental matrix $M(t)$ of the solution on $[0,2\pi]$ can be obtained (by the standard Piquard iterations) as the sum $$M(t)=\begin{pmatrix} 1 & 0 \cr 0 & 1 \end{pmatrix}+ \begin{pmatrix} 0 & \psi_1(t) \cr -\psi_1(-t) & 0 \end{pmatrix}$$ $$+ \begin{pmatrix} -\psi_2(t) & 0 \cr 0 & -\psi_2(-t) \end{pmatrix}+ \begin{pmatrix} 0 & -\psi_3(t) \cr \psi_3(-t) & 0 \end{pmatrix}+ O(C^4)$$ Where $\psi_0(t)=1$ and $\psi_{k+1}(t)=\int_0^t\psi(s)\psi_k(-s)ds$.

You can write a long series but I want to convince you that 4 first terms are enough for your problem. We can find the first two $\psi$'s:

$$\psi_1=\frac 1i[(ce^{it}-\bar c{e^{-it}})-(c-\bar c)$$

$$\psi_2=-\frac{c-\bar c}{i}\psi_1+\frac 1i(c^2-\bar c^2)t+\frac{|c|^2}2(e^{2it}+e^{-2it})-|c|^2$$

and the linear term in $\psi_3$, which is $$(c^2-\bar c^2)t[(c-\bar c)+ce^{it}-\bar ce^{-it}]$$

Plugging in $t=2\pi$, we see that the monodromy matrix is $\begin{pmatrix} 1-2\pi v & 4\pi vs \cr 4\pi vs & 1+2\pi v \end{pmatrix}+\begin{pmatrix} O(C^4) & O(C^5) \cr O(C^5) & O(C^4)\end{pmatrix}$ with $v=2\Im (c^2)=\frac 12 C^2\sin 2\Phi= \frac 12 \alpha^{-2}AB\sin(\Phi_B-\Phi_A)$, $s=c-\bar c=C\sin\Phi=\alpha^{-1}\sqrt{AB}\sin\frac{\Phi_B-\Phi_A}2$.

Now, the life is easy: the growth/decay part is essentially given by the matrix $$G(t)=\begin{pmatrix} e^{-vt} & 0 \cr 0 & e^{vt} \end{pmatrix}$$ the rotation part is essentially given by $$T=\begin{pmatrix} 1 & -s \cr s & 1 \end{pmatrix}$$ and the oscillation is essentially given by $$H(t)=\begin{pmatrix} 1-\widetilde\psi_2(t) & \psi_1(t) \cr -\psi_1(-t) & 1-\widetilde\psi_2(-t) \end{pmatrix}$$ where $\widetilde \psi_2$ is $\psi_2$ with the term $\frac 1i(c^2-\bar c^2)t$ removed. So, $M(t)\approx H(t)T^{-1}G(t)T$ (well, $T^{-1}GT$ and $H$ don't commute but the commutator effect is of size $C^3$).

This should work just fine letting you to see just enough in your range.

P.S. My original answer had an error in that I neglected the rotation of eigenvectors, but now it should be fine even for fairly large $C$ like $0.1$. Check agaist your numerics and see if it works well enough for you.

-
Hi, I hope you read still this, but with a project proposal, holidays and last but not least reading into the background of your answer, it essentially became ... now :) I have two questions left: 1. I do understand now why the monodromy matrix is self-adjoint, but what do you derive from that, or why is it important? 2. How do you come from the monodromy matrix to G/T/H and thus the approximated M? I think I am just blind here. –  elcron Aug 25 '11 at 14:16
1) It is not terribly important but quite convenient because the eigenvectors are orthogonal, etc. 2) The approximation is "exact" at $2\pi k$ (up to high order terms). In between it would be just as good if the matrices commuted. They don't, but, if you look at the commutator, it is of higher order. You see, essentially I just take the $t$-th power of the monodromy matrix ($T^{-1}G(t)T$) and correct for the periodic "phase factor" $H(t)$. They interact in some fancy way due to non-commutativity, but I just ignore it in the low order approximation. –  fedja Aug 26 '11 at 1:33
To 1) But you don't use that anywhere, do you? To 2) I think might be missing some basic theorems here again (as with the self-adjoined part :)). H contains the first 3 parts of the result obtained by iteration. As the 4th part is not periodic anymore, you don't include it in the oscillation part. Correct? Now, you still include it into the monodromy matrix, that governs the grows (exact to the iteration at $2\pi k$). Correct? What I still don't understand: How do you go from the monodromy matrix to T and how is $T^{-1}G(t)T$ the t-th power? What books do I have to look into for that? –  elcron Aug 26 '11 at 7:04
I just find the eigenvectors and the eigenvalues (up to high order terms) and do the usual diagonalization. –  fedja Aug 28 '11 at 10:29
Ok, I am starting to understand. I can reproduce the monodromy matrix (s should read $(c-\bar c)/i$, right?) and a 2nd order Taylor series of the eigenvalues look like G. Is there a more exact way to see that? I still see no light in getting that the eigenvalues are (1,s),(-s,1), but if you say so there should be a way :) Currently I am drowning in complex conjugates and stuff. Why do you stop at $\psi_2$ for H(t)? Also I would still apprciate a book recommendation. I looked several about ODEs (Taylor, Arnold,Farlow,...) but they hardly even mention monodromy. –  elcron Aug 29 '11 at 9:30
show 2 more comments

Set $\alpha z=w$, you get the new system $$dy/dw={1\over \alpha}\pmatrix{0&B\cos(w+\Phi_b)\cr A\cos(w+\Phi_a)&0}y.$$ Since you are interested in a case where $\alpha$ is large and $w$ is of moderate size, you can try an expansion of the solution in powers of $1/\alpha$.

-
add comment

Another thing to consider: introduce new independent variable $\tau=(1/\alpha)sin(\alpha z+\Phi_a)$. Also note that we have $$cos(\alpha z+\Phi_b)=cos(x+\delta)$$ and $$cos(x+\delta)=cos(x)cos(\delta)-sin(x)sin(\delta).$$ where $x=\alpha z+\Phi_a$ and $\delta=\Phi_b-\Phi_a$.

Then the first equation of your system becomes $$dy_1/d\tau=B \left(cos(\delta)-\displaystyle\frac{\alpha\tau}{\sqrt{1-\alpha^2\tau^2}} sin(\delta)\right) y_2,$$ while the second one is simply $$dy_2/d\tau=A y_1.$$ This new system might be somewhat easier to investigate, be it analytically or numerically.

-
add comment