Coupled linear time-varying system of differential equations. Analytical solution?

Question

How to solve the system analytically in the interval {to,t}

x''[t]==-a0*(a1 - b*z'[t])*Cos[w*t]

z''[t]=-a0*b*x'[t]*Cos[w*t]

This gives an error:

eqn1 = {x''[t] == -a0*(a1 - b*z'[t])*Cos[w*t], x[t0] == 0, x'[t0] == 0};
eqn2 = {z''[t] == -a0*b*x'[t]*Cos[w*t], z[t0] == 0, z'[t0] == 0};

soln = DSolve[{eqn1, eqn2}, {x,z}, t][[1]]

This is free electron in electromagnetic field. $m_e\begin{bmatrix}x''\\y''\\z''\end{bmatrix}=-e\begin{bmatrix}E\\0\\0\end{bmatrix}Exp(iwt)-e\begin{bmatrix}x'\\y'\\z'\end{bmatrix}$x$\begin{bmatrix}0\\B\\0\end{bmatrix}Exp(iwt)$=$-e[\begin{bmatrix}E\\0\\0\end{bmatrix}+\begin{bmatrix}-z'B\\0\\x'B\end{bmatrix}]Exp(iwt)$

If the fields are complex $Exp(iwt)$ mathematica finds some solution to the reduced equation of 3rd order for $z$. But the answer is also complex. Is it possible to get real Z for cos(wt) field out of it?

eqns = {z'''[t] == a*Exp[I*t*2*w]*(1 - z'[t]*g) + z''[t]*(I*w), 
   z[t0] == 0, z'[t0] == 0};

soln = Real[DSolve[eqns, z, t][[1]]]

$a=e^2*E^2/(m_e^2c)$; $ g=m_e/(e*c)$; $w=2\pi/\lambda$, where $B=E/c$

Let say: $t0$ inside $(0,pi/w)$ , $t$ is $(0,2\pi/w)$, $w=2π/λ$

Let say $E=10^9V/m$, $B=E/c$, electron $m= 9×10^{-31}kg$, $e = -1.602×10^{−19}C$;

Answer 1

Here is an analytical solution. The solution is based on series expansion of BesselJ functions of first kind.

It turned out the the solution to

\begin{align*} x^{\prime\prime}\left( t\right) & =-a_{0}\left( a_{1}-bz^{\prime}\left( t\right) \right) \cos\left( \omega t\right) \\ z^{\prime\prime}\left( t\right) & =-a_{0}bx^{\prime}\left( t\right) \cos\left( \omega t\right) \end{align*}

With initial conditions $x\left( 0\right) =0,x^{\prime}\left( 0\right) =0,z\left( 0\right) =0,z^{\prime}\left( 0\right) =0$ is

\begin{align*} x\left( t\right) & =\frac{-a_{1}}{b}\int_{0}^{t}\sin\left( \frac{a_{0}% b}{\omega}\sin\left( \omega x\right) \right) dx\\ z\left( t\right) & =\frac{1}{b}\int_{0}^{t}a_{1}\left( 1-\cos\left( \frac{a_{0}b}{\omega}\sin\left( \omega x\right) \right) \right) dx \end{align*}

For $\omega\neq0$. (Thanks to Maple 2016 help). However, how to evaluate the above? It turned out that, with the help of looking up BesselJ identities and googling around that

$$ \sin\left( \frac{a_{0}b}{\omega}\sin\left( \omega x\right) \right) =2\sum_{m=0}^{\infty}J\left( 2m+1,\frac{a_{0}b}{\omega}\right) \sin\left( \left( 2m+1\right) \omega x\right) $$

And

$$ \cos\left( \frac{a_{0}b}{\omega}\sin\left( \omega x\right) \right) =J\left( 0,\frac{a_{0}b}{\omega}\right) +2\sum_{m=1}^{\infty}J\left( 2m,\frac{a_{0}b}{\omega}\right) \cos\left( 2m\omega x\right) $$

And these can be analytically integrated by Mathematica! It also turned out, that taking the limit to small value is enough to get excellent result. For example, $m=10$. Using the above, here is the full analytical solution for $m=10$ and comparing it to the numerical solution also

\begin{align*} x\left( t\right) & =\frac{-a_{1}}{b}\int_{0}^{t}2\sum_{m=0}^{\infty }J\left( 2m+1,\frac{a_{0}b}{\omega}\right) \sin\left( \left( 2m+1\right) \omega x\right) dx\\ z\left( t\right) & =\frac{1}{b}\int_{0}^{t}a_{1}\left( 1-\left[ J\left( 0,\frac{a_{0}b}{\omega}\right) +2\sum_{m=1}^{\infty}J\left( 2m,\frac{a_{0}% b}{\omega}\right) \cos\left( 2m\omega x\right) \right] \right) dx \end{align*}

The rest now is show using Mathematica

Clear[a0, b, a1, w, t0, tf];
f[n_, t_] := 2 Sum[BesselJ[2 m + 1, a0 b/w] Sin[(2 m + 1)  w t], {m, 0, n}];
solX = (-a1/b) Integrate[f[10, x], {x, 0, t}]

params={a0->2,b->3, a1->3,w->2};
t0=0;tf=1;
pX=Plot[solX/.params,{t,t0,tf},ImageSize->300,PlotLabel->"Analytical solution of x(t)"]

For z(t)

g[n_,t_]:=BesselJ[0,a0 b/w]+2 Sum[BesselJ[2m,a0 b/w] Cos[2m  w t],{m,1,n}]
solZ=(1/b)Integrate[a1(1-g[10,x]),{x,0,t}]
t0=0;tf=1;
pZ=Plot[solZ/.params,{t,t0,tf},ImageSize->300,PlotLabel->"Analytical solution of z(t)"]

Compare to Numerical

Clear[x,z,t,w,a0,b,t0,tf,m]
a0=2;b=3; a1=3;w=2;t0=0;tf=1;
x''[t]==-a0*(a1-b*z'[t])*Cos[w*t]
z''[t]==-a0*b*x'[t]*Cos[w*t]

eqn1={x''[t]==-a0*(a1-b*z'[t])*Cos[w*t],x[t0]==0,x'[t0]==0};
eqn2={z''[t]==-a0*b*x'[t]*Cos[w*t],z[t0]==0,z'[t0]==0};

soln=NDSolve[{eqn1,eqn2},{x,z},{t,t0,tf}]

pXN=Plot[x[t]/.soln,{t,0,1},PlotLabel->"x(t) Numerical",ImageSize->300]
pzN=Plot[z[t]/.soln,{t,0,1},PlotLabel->"z(t) Numerical",ImageSize->300]

Putting them on same plot, shows they are the same

Grid[{Show[pX,pXN],Show[pZ,pzN]}]

P.S. Maple solution I used as starting point is below

restart;
eq1:=diff(x(t),t$2)= -a0*(a1-b*diff(z(t),t))*cos(w*t);
ic1:=x(0)=0,D(x)(0)=0;
eq2:=diff(z(t),t$2)= -a0*b*diff(x(t),t)*cos(w*t);
ic2:=z(0)=0,D(z)(0)=0;
dsolve({eq1,eq2,ic1,ic2},{x(t),z(t)}) assuming(w>0);

Answer 2

NDSolve is probably a better option:

ClearAll["Global`*"]

e = 1.60217662*^-19;(*coulombs*)
m = 3.343583719*^-27;(*D+kg*)

eF = 100;(*V/m*)
bF = 0.01;(*Tesla*)

w = 1;

tmax = 1*^-4;

nds = NDSolve[{
    x''[t] == -e/m*Cos[w*t]*(eF - bF*z'[t]),
    z''[t] == -e/m*bF*Cos[w*t]*x'[t],
    x[0] == 0, z[0] == 0, x'[0] == 0, z'[0] == 0},
   {x, z}, {t, 0, tmax}];

ParametricPlot[Evaluate[{z[t], x[t]} /. nds], {t, 0, tmax}]

Take a look at this post as well: Using NDSolve to find particle trajectory

Update

I don't think there is a closed form solution available from DSolve:

ClearAll["Global`*"]
opts = Options[Solve];
SetOptions[Solve, Method -> Reduce];
DSolve[{
   x''[t] == -e/m*Cos[w*t]*(eF - bF*z'[t]),
   z''[t] == -e/m*bF*Cos[w*t]*x'[t],
   x[0] == 0, z[0] == 0, x'[0] == 0, z'[0] == 0},
  {x[t], z[t]}, t] // FullSimplify
SetOptions[Solve, opts];

This returns "un-evaluated". For challenging ODE analytic solutions take a look at Why does DSolve return two solutions for my ODE? and Accessing Reduce from DSolve

Update-2

You updated the question with essentially the same equations so my previous comments/answer still applies. For completeness here is the updated equation that returns un-evaluated.

eqn1 = {x''[t] == -a0*(a1 - b*z'[t])*Cos[w*t], x[0] == 0, x'[0] == 0};
eqn2 = {z''[t] == -a0*b*x'[t]*Cos[w*t], z[0] == 0, z'[0] == 0};

opts = Options[Solve];
SetOptions[Solve, Method -> Reduce];
DSolve[Join[eqn1, eqn2], {x[t], z[t]}, t] // FullSimplify 
SetOptions[Solve, opts];

Answer 3

It is possible to obtain an analytic solution, at least for the derivatives.

eqn1 = {x''[t] == -a0*(a1 - b*z'[t])*Cos[w*t], x[t0] == 0, x'[t0] == 0};
eqn2 = {z''[t] == -a0*b*x'[t]*Cos[w*t], z[t0] == 0, z'[t0] == 0};

In principle, since the equations are linear in x and z it should be straight forward to achieve any required initial conditions by adding a linear function of time, so I'll ignore the initial conditions.

Make the substitution

subst = {X[t] == x'[t], Z[t] == z'[t]}

converting from an apparently second order differential equations to first order equations.

neweqn = {First[eqn1], First[eqn2]} /. {x -> X, z -> Z} /. Derivative[n_] -> Derivative[n - 1]

This has an analytic solution

sol1 = DSolve[neweqn, {X, Z}, t]
(* {{X -> Function[{t}, C[1]*Cos[(a0*b*Sin[t*w])/w] + C[2]*Sin[(a0*b*Sin[t*w])/w]], 
Z -> Function[{t}, a1/b + C[2]*Cos[(a0*b*Sin[t*w])/w] -   C[1]*Sin[(a0*b*Sin[t*w])/w]]}} *)

whose form, I imagine, could be enlightening.

We can proceed to a solution in the original functions

sol2 = DSolve[subst /. sol1, {x, z}, t]

but this gives an integral that Mathematica cannot simplify. We can verify that this is a solution of the original equations

Flatten[First /@ {eqn1, eqn2} /. sol2] // Simplify
(* {True, True} *)

Answer 4

The question already contains the analytic solution for z with complex-exponential driving, so I think the only remaining question is how to get the solution for cosinusoidal driving.

First I re-derive the third-order equation quoted in the question from the initial set of two equations (with unit chosen such that $e=1$ and $E=a$):

eqx = x''[t] == -(a - b*z'[t])*Exp[I w*t];
eqz = z''[t] == -b*x'[t]*Exp[I w*t];
eq = Eliminate[D[eqz, t] && eqz && eqx, {x'[t], x''[t]}]

$$i \left(a b e^{2 i t w}+b^2 \left(-e^{2 i t w}\right) z'(t)-z^{(3)}(t)\right)=w z''(t)$$

Then the exact analytical solutions can be obtained if the initial conditions for $z(t_0)$, $z'(t_0)$ and $x'(t_0)$ are specified (the latter was omitted in the question). I denote the initial velocity $x'(t_0)$ by vx0. The solution for complex driving is defined as a function of the parameters and time:

Clear[s];
s[w_, t0_, a_, b_, vx0_] = 
  z /. First@
    DSolve[{eq, z[t0] == 0, z'[t0] == 0, 
      eqz /. t -> t0 /. x'[t0] -> vx0}, z, t];

Clear[sCos];
sCos[w_, t0_, a_, b_, vx0_][t_] := (s[w, t0, a, b, vx0][t] + 
                                     Conjugate[s[w, t0, a, b, vx0][t]])/2

The second definition is the solution to the same problem with cosinusoidal driving. All I did was to superimpose the solution to the exponential driving problem with its complex conjugate, and divided by 2. This works because the original system of equations for cosinusoidal driving is linear and has entirely real coefficients. I use the Conjugate instead of simply reversing the sign of the frequency w in the superposition, because the DSolve solution may contain branch cuts which could make the result of $\omega\to-\omega$ not equal to the complex conjugate of the function. Here are some examples of how to evaluate the solutions (s for exponential and sCos for cosinusoidal driving):

s[1., 0, 1, .1, -1.1][10]

(* ==> 125.918 - 0.052946 I *)

sCos[1., 0, 1, .1, -1.1][10]

(* ==> 125.918 + 0. I *)

To avoid having Conjugate appear in the solution, you could also define

Clear[s2];
s2[w_, t0_, a_, b_, vx0_] = 
  z /. First@
    DSolve[{eq /. w -> -w, z[t0] == 0, z'[t0] == 0, 
      eqz /. t -> t0 /. x'[t0] -> vx0}, z, t];

Clear[sCos];
sCos[w_, t0_, a_, b_, vx0_][t_] := 
         (s[w, t0, a, b, vx0][t] + s2[w, t0, a, b, vx0][t])/2