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

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$;

• There appears to be a space or * missing between E and e. Fixing that, and removing everything outside your DSolve[...] shows that DSolve can't find a solution and returns the original which is not of the form {x->stuff,y->stuff} and this is what gives you your error message when you try to use it with /.
– Bill
Jul 8, 2016 at 22:37
• First consider only the DSolve term. Mathematica can't do that. It may be able to solve it if you do a substitution x'[t]->u[t] and z'[t]->v[t]' Jul 8, 2016 at 22:39
• Thanks. How to make it solve then DSolve[x''[t] == -Cos[w*t]*eo, x[y] == 0, x'[y] == 0, x[t], t] , y is particular t Jul 8, 2016 at 22:48
• What is you anticipated time scale values? And values for the other constants? Jul 9, 2016 at 16:31
• $t0 = (0,pi/w)$, $w=2\pi/\lambda$, t is (0,2\pi/w), Let say E=10^9V/m, B=E/c, electron m= 9e-31, e = -1.602×e^−19; Jul 9, 2016 at 18:11

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);


• Very interesting +1. I'm having trouble matching the results to the numeric solution when the variable values are close to "real world variables" - params = {a0 -> e/m, b -> bF, a1 -> eF, w -> w} and tf=tmax. I think the issue may be with tf being small. Jul 9, 2016 at 13:12
• It seems the analytical equations are having trouble matching the cycloid motion (not just because of variable scale). Jul 9, 2016 at 13:26
• $n$ might need to be impractically large Jul 9, 2016 at 13:36
• @Young are right, for the very small time scale you are trying it on, you'll need a very large sum of the BesselJ series. I tried up to 50,000 terms and started to get close to the numerical, but it got very slow by then. As always, sometimes series solution might need very large number of terms to converge. Better stick to the numerical solution for such small time scale. For large time scale, you can get away with small number of terms it seems. Jul 9, 2016 at 15:13
• Mathematica is solving the case where the field is comlex instead of cos. It is same equation still if I take the real part at the end? eqns = {z'''[t] == a0*Exp[It*2*w]*(1 - z'[t]*a1) + z''[t]*(Iw), z[t0] == 0, z'[to] == 0}; soln = Real[DSolve[eqns, z, t][[1]]] Jul 9, 2016 at 17:29

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];


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} *)


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

• Thanks how to let mathematica know that the constants are real ? Jul 9, 2016 at 22:08
• You could do Assuming[{w,t0,a,b,vx0}\[Element]Reals,Simplify@ComplexExpand[sCos[w,t0,a,b,vx0][t]]] but it doesn't get rid of all the imaginary units. The CoshIntegral` etc. appearing in the solutions appear to be too hard to turn into purely real functions. No matter what you do, I think the analytic solution will always look unwieldy.
– Jens
Jul 9, 2016 at 22:19