# A hint for an exact solution of a system of differential equations

I am looking for a curve in two-dimensional space, so $c(t)=(x(t),y(t))$, that satisfies a system of differential equations which could not be solved with Mathematica $$x''(t) = (k - c \; y(t))\; y'(t) \\ y''(t) = -(k - c \; y(t))\; x'(t)$$ The differential equations enforce that $c''(t)$ is perpendicular to $c'(t)$, so it is clear that $c'(t)$ always has constant length and it may be assumed this length is $1$.

The constants $k$ and $c$ are both positive. I am interested especially in the case where $k$ is about $1/4$ and $c$ is about $4/1000$, but I think it can be solved with arbitrary positive $k$ and $c$.

The differential equations can be solved quite good numerically. And I know from numerical computations that a prolate cycloid has nearly the same orbit (path) as the seeked curve. Unfortunately the solution is only nearly a prolate cycloid. One deficiency of a prolate cycloid is that it doesn't have constant velocity $|c'(t)|$.

I would be very glad, if someone could give me a hint to solve these equations.

• Can you please post your code? That would provide folks more info on what may have gone wrong. Feb 23, 2017 at 13:12
• Also, do not use uppercase names for your variables, especially for single letter ones, as they may conflict with built in symbols. In your case, both K and C do in fact have built in meaning; it is safer if you replace them with lowercase. Feb 23, 2017 at 13:16
• What are the four conditions? You mean c(t) is a phase-space trajectory?
– zhk
Feb 23, 2017 at 13:17

Symbolic solution

Even though Mathematica cannot solve straightforwardly the original system we can get the symbolic solutions with a little bit smarter approach. Let's observe that we can get rid of x'[t] and x''[t] from the system by changing variables y -> z[t] == k - c y[t], thus z'[t] == -c y'[t] and z''[t] == -c y''[t] therefore x''[t] == (k - c y[t]) y'[t] implies that x''[t] == -1/c z[t] z'[t], integrating it we get c x'[t] + 1/2 z[t]^2 + d == 0, where d is a constant of integration. Inserting this into the second equation yields

    z''[t] + 1/2 z[t]^3 + d z[t] == 0


This equation can be solved in terms of Jacobi elliptic functions:

sols = z[t] /. DSolve[{z''[t] + 1/2 z[t]^3 + d z[t] == 0}, {z[t]}, {t}]


Recalling that z[t] == k - c y[t] we can easily figure out periodic-like behaviour of y[t] found in numerical calculations. In order to get x[t] we need inserting y[t] into the original system, then integrating with respect to t. It seems that we should work with appropriate initial conditions y[0] and y'[0] and choose adequate branches of solutions. Then various special cases could clarify the overall behaviour of the system, let's consider a special case.

 sd041 = FullSimplify[sols /. {d -> 0, C[1] -> 4, C[2] -> 1}, t > 0]

{-2 JacobiSN[1 + t, -1], 2 JacobiSN[1 + t, -1]}


Taking the second equation of the original system we can get x'[t] (assuming c == 1 and k == 1):

FullSimplify[ D[ sd041[[1]], {t, 2}]/(1 - sd041[[1]]), t > 0]

(4 JacobiSN[1 + t, -1]^3)/(1 + 2 JacobiSN[1 + t, -1])


and x[t] is (putting the constant of integration equal to zero):

xx = Integrate[(4 JacobiSN[1 + t, -1]^3)/(1 + 2 JacobiSN[1 + t, -1]), t]


We plot only the real part of xx as x[t], (this is because of another issue, see e.g. Why does Integrate declare a convergent integral divergent? )

Plot[{ Re @ xx, 1 - sd041[[1]]}, {t, 0, 20}]


and the solution in the phase space is

ParametricPlot[{ Re @ xx, 1 - sd041[[1]]}, {t, 0, 20}]


At this point it is plausible to exploit numerical capabilities of the system

Numerical solution

We can exploit NDSolve and including arbitrary initial conditions, we can compare various cases:

With[{c = 1/4, k = 4/1000},
ds = NDSolve[{x''[t] == (k - c y[t]) y'[t],
y''[t] == -(k - c y[t]) x'[t],
x[0] == 1, y[0] == 0, x'[0] == 2, y'[0] == -1},
{x[t], y[t]}, {t, 0, 150}]];
{X[t_], Y[t_]} = {x[t], y[t]} /. First @ ds;


Now we plot the solution in the range 0 < t < 50:

Plot[{X[t], Y[t]}, {t, 0, 50}]


as well as its derivative

Plot[{X'[t], Y'[t]}, {t, 0, 50}]


Another interesting feature can be observed with the parametric plot of the solution:

ParametricPlot[{X[t], Y[t]}, {t, 0, 60}]


and it is expected to compare the first and the second derivatives of the solution, we can do it with the animated parametric plot:

tab =
Table[ ParametricPlot[{{X'[t], Y'[t]}, {X''[t], Y''[t]}}, {t, 0, v},
PlotLegends -> Placed[Style[Row[{"t = ", NumberForm[N@v, {3, 1}]}],
Bold, 20], {Left, Top}],
PlotRange -> {{-3.5, 3.5}, {-3.5, 3.5}}], {v, 1, 16, 1/6}];


ListAnimate[ tab, Paneled -> False]

Edit

The OP expected solutions with reversed values (c == 4/1000, k == 1/4) Then the system behaves in a different way:

With[{c = 4/1000, k = 1/4},
ds = NDSolve[{x''[t] == (k - c y[t]) y'[t],
y''[t] == -(k - c y[t]) x'[t],
x[0] == -16, y[0] == 4,
x'[0] == 12, y'[0] == -4}, {x[t], y[t]}, {t, 0, 130}]];
{X[t_], Y[t_]} = {x[t], y[t]} /. First@ds;


Now the phase space solution is more similar to cycloid-like

ParametricPlot[{X[t], Y[t]}, {t, 0, 100}]


 ParametricPlot[{{X'[t], Y'[t]},
{X''[t], Y''[t]}}, {t, 0, 100},
PlotLegends -> "Expressions"]


Playing with different initial conditions we can get various patterns of the behaviour, more or less similar to the previous case with {c = 1/4, k = 4/1000}. Nonetheless from programming point of view the task of examining the system is the same.

• In OP's question K=1/4 and C=4/1000, why you swap it? BTW, thanks for the nice answer.
– zhk
Feb 23, 2017 at 15:28
• @MMM Thanks for pointing out that, indeed now the solution looks a bit different, nevertheless I suppose that imposing an appropriate initial data, we could get more similar solution. The OP has not mentioned what initial condition are needed. Feb 23, 2017 at 16:33
• Thank you Artes for your nice numerical solutions. The patterns that come out are really beautiful. But I mentioned that I already know how to get a numercal solution. What I am looking for is an analytical solution, that means to have a formula which solves the equation. Feb 23, 2017 at 16:55
• Please provide more detail on obtaining z''[t] + 1/2 z[t]^3 + d z[t] == 0. Thanks. Feb 23, 2017 at 20:11
• @Miene There is no claim in the OP that you want to find a symbolic solution. Edit your question, and write unambigously what you've expected. Feb 24, 2017 at 11:05

Analytical solution with Maple

Maple is able to produce an analytical solution with out the use of the initial conditions,

restart;
Eq1:= diff(x(t),t$2)=(k-c*y(t))*diff(y(t),t$1):
Eq2:= diff(y(t),t$2)=-(k-c*y(t))*diff(x(t),t$1):
dsolve({Eq1,Eq2});


But, if when we specify only two initial conditions ,

sol1:=dsolve({Eq1,Eq2,x(0)=a1,y(0)=b1});


odetest(sol1,{Eq1,Eq2,x(0)=a1,y(0)=b1})


{0}

If we can take another combination of the conditions,

sol2:=dsolve({Eq1,Eq2,x(0)=a1,D(x)(0)=c1});


odetest(sol2,{Eq1,Eq2,x(0)=a1,D(x)(0)=c1})


{0}

With these combinations y(0)=b1,D(y)(0)=d1 and D(x)(0)=c1,D(y)(0)=d1 maple produces no output.

Mathematica's DSolve is unable to solve the system in question analytically.

Numerical solution with Mathematica

C1 = 4/1000; K1 = 1/4;

sol1[x0_?NumericQ] :=
First@NDSolve[{x''[t] == y'[t]*(K1 - C1*y[t]),
y''[t] == -x'[t]*(K1 - C1*y[t]), x[0] == x0, x'[0] == x0,
y[0] == x0, y'[0] == x0}, {x, y}, {t, -10, 10}];

ParametricPlot[
Evaluate[{x[t], y[t]} /. sol1[#] & /@ Range[-10, 10, 1]], {t, -10,
10}, Frame -> True, PlotRange -> All]


First I appreciate the revision of my initial question done by Artes. It is much clearer now. Further I thank all who gave me hints for analytical solutions. In the end I had to accept that there is no really simple solution --- as I posted my question, I had the hope there might be a solution which I just can see even it is not to intricate.

After all I had to go back to a numerical solution and I'll describe my proceeding as it is strait forward and can even be carried out on Excel:

I start with a initial condition c(0) and c'(0) [c'(0) with length 1] and the differential equation gives me c''(0) as

x′′(0)=(k−cy(0))y′(0)
y′′(0)=−(k−cy(0))x′(0).


Now for a small t I get

x'(t)=x'(0)+t*x''(0)
y'(t)=y'(0)+t*y''(0).


I normalize c'(t) to length 1 and take the mean value mv(t) = (c'(0)+c'(t))/2, which I normalize likewise. Then I count

x(t)=x(0)+t*mvx(t)
y(t)=y(0)+t*mvy(t),


where mvx is the x-part of mv an mvy is the y-part of mv. Then I do the same procedure with starting point c(t) and c'(t) to get values for 2*t and so on ...

In this way I get c(n*t) for n=0, 1, 2, 3, ... If t is small enough this will be a fine approximation to the real solution.