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.

  • 4
    $\begingroup$ Can you please post your code? That would provide folks more info on what may have gone wrong. $\endgroup$
    – dearN
    Feb 23, 2017 at 13:12
  • 3
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    Feb 23, 2017 at 13:16
  • 2
    $\begingroup$ What are the four conditions? You mean c(t) is a phase-space trajectory? $\endgroup$
    – zhk
    Feb 23, 2017 at 13:17

3 Answers 3


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

enter image description here

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

enter image description here

and the solution in the phase space is

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

enter image description here

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

enter image description here

as well as its derivative

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

enter image description here

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

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

enter image description here

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]

enter image description here


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

enter image description here

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

enter image description here

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.

  • $\begingroup$ In OP's question K=1/4 and C=4/1000, why you swap it? BTW, thanks for the nice answer. $\endgroup$
    – zhk
    Feb 23, 2017 at 15:28
  • $\begingroup$ @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. $\endgroup$
    – Artes
    Feb 23, 2017 at 16:33
  • $\begingroup$ 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. $\endgroup$
    – Miene
    Feb 23, 2017 at 16:55
  • 1
    $\begingroup$ Please provide more detail on obtaining z''[t] + 1/2 z[t]^3 + d z[t] == 0. Thanks. $\endgroup$
    – bbgodfrey
    Feb 23, 2017 at 20:11
  • 1
    $\begingroup$ @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. $\endgroup$
    – Artes
    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,

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

enter image description here

But, if when we specify only two initial conditions ,


enter image description here



If we can take another combination of the conditions,


enter image description here



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

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

enter image description here


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


Now for a small t I get


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


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.