@Rouben Rostamian  

But in this particular case, the classical ode has a C^1 solution which is C^2 except a finite set of points (in a bounded interval).
This solution obviously satisfies the differential inclusion.

@Preben Alsholm 

The formula

sol:=(A-n*C)*cos(omega*t)+(-1)^n*C/2;

is valid only for t < tmax = A*Pi*sqrt(k/m)/(2*mu*g)

because after that   A - nC  becomes <0.
That is why the graph obtained by dsolve is different from sol for t>tmax.

Another interesting thing. sol is only an approximation! A rather good one.
It is possible to compute the exact solution but this will be of the form

A(t)*cos(omega*t + alpha(t)) + S(t)  where A(t), alpha(t) , S(t) are locally (piecewise) constant.

I like this problem!   :-)

Edit. I was intrigued and I have computed by hand the solution. sol is correct (i.e. alpha(t)=0), it is dsolve which makes the difference (for t<tmax of course).

@Preben Alsholm 

The following graph (yours) includes (-1)^n in red. It is obvious now that the ode refers to signum(x') not signum(x'');

@tomleslie 

It depends on how the ode is understood. A natural assumption is that x''() exists and is nonzero except on a discrete set A.
In this case x() cannot be continuous.
Otherwise any piecewise linear function is a solution because x''() = 0 except on a discrete set, so we don't have any condition for x(). So, it can even oscillate:

@Preben Alsholm 

The ode is given but the problem is that x() cannot be continuous.
So, the desired shape of the solution is out of the question.

This direct consequence of Green's formula is used:

     ∫∫_D f(x,y) dx dy = ∫_C Q(x,y) dy,  where Q := ∫ f(x,y) dx

It remains to compute the curvilinear integral.
We must find a positively oriented parametrization (by hand!):

x=X(t), y=Y(t), t ∈ [a,b] of the (closed) curve C.

Then the curvilinear integral reduces to the usual Riemann integral

       ∫_a^b Q(X(t),Y(t))·Y^′(t) dt

The boring (sometimes difficult) part is to write the parametrization (usually piecewise).

@Kitonum 

I have used a compiled procedure.
BTW, I have corrected the comment which was wrong.
X(10^9)  needed several minutes but X(100000) only 65ms.

@rlopez 

It is solve which returns infinitely many solutions. Yes, if only the principal branch of ln is to be considered (as Maple does) then there are only two solutions, because the equation is equivalent to  ln(z) = k*Pi + k*Pi*I  (k in Z), which is possible only for k=0 and k=1 (argument(z) = k*Pi must be in (-Pi,Pi]).

# X(10^7) = 29529388;
# X(10^8) = 213554980;
# X(10^9) = 92599500093871447;

(Edited)

@Earl 

You may fine-tune the dsolve parameters. E.g. adding method=bvp[midrich]
is enough to work for a smaller speed.

@Christopher2222 

The following list has a simple and clear pattern. But probably it is (almost) impossible to find (guess) the next number.

1, 4, 10, 19, 31, 41, 49, 59, 78, 109, 143, 166, 173, 178, 204, 259, 328, 378, 388, 377, 388, 457, 570, 668, 701, 672, 647, 703, 850, 1018, 1106, 1078, 1004, 1009, 1161, 1401, 1582, 1596, 1483, 1403, 1509, 1797, 2096, 2210, 2097, 1922, 1923, 2201, 2613, 2882, 2841, 2594, 2445, 2634, 3111, 3565, 3681, 3433, 3124, 3137, 3592, 4215, 4566, 4421, 3996, 3773, 4087, 4806, 5429, 5511, 5072, 4605, 4659, 5346, 6215, 6628, 6325, 5678, 5390, 5881, 6893, 7688, 7687, 7002, 6362, 6498, 7477, 8618, 9061, 8538, 7632, 7301, 8031, 9382, 10338, 10197, 9210, 8396, 8665, 9999, 11426, 11856, 11045, 9853, 9513, 10551, 12280, 13375, 13027, 11686, 10707, 11174, 12923, 14642, 14999, 13834, 12336, 12037, 13455, 15596, 16793, 16161, 14420, 13300, 14040, 16263, 18263, 18479, 16889, 15080, 14883, 16758, 19334, 20580, 19585, 17405, 16183, 17277, 20028, 22286, 22280, 20200, 18085, 18066, 20474, 23497, 24726, 23282, 20635, 19365, 20903, 24226, 26704, 26386, 23754, 21355, 21601, 24618, 28084, 29217, 27237, 24106, 22857, 24933, 28866, 31508, 30781, 27541, 24897, 25503, 29202, 33092, 34037, 31436, 27818, 26673, 29384, 33950, 36687, 35446, 31554, 28718, 29792, 34237, 38517, 39168, 35864, 31772, 30830, 34272, 39482, 42226, 40363, 35786, 32829, 34486, 39733, 44351, 44594, 40508, 35973, 35344, 39613, 45461, 48110, 45517, 40234, 37245, 39603, 45697, 50583, 50293, 45358, 40427, 40236, 45421, 51883, 54320, 50888, 44898, 41982, 45163, 52133, 57199, 56246, 50405, 45146, 45524, 51708, 58743, 60837, 56462, 49778, 47057, 51183, 59046, 64183, 62432, 55641, 50141, 51232, 58486, 66031, 67640, 62225, 54881, 52492, 57683, 66433, 71518, 68832, 61063, 55428, 57380, 65764, 73735, 74705, 68162, 60214, 58308, 64677, 74291, 79182, 75425, 66669, 61026, 63990, 73547, 81840, 82011, 74263, 65788, 64528, 72182, 82615, 87152, 82194, 72461, 66956, 71084, 81838, 90328, 89533, 80520, 71618, 71178, 80211, 91393, 95405, 89120, 78445, 73240, 78683, 90639, 99178, 97248, 86927, 77721, 78281, 88773

@Rouben Rostamian  

Yes, I made a mistake, unfortunately I did not check the statement about uniquness. y(1-x) is a path joining the points (0,0) and  (1,0) but not an admissible one (the boat will not have a constant speed here).
Also, because -Pi/2 < alpha, psi < Pi/2, it is not difficult to see that psi-alpha = Pi - arcsin(k)   cannot occur.
Sorry, my excuse is that I'm on holyday :-)

The optimal path is not unique in general: if y(x) is an optimal path, then y(1-x) is optimal too because the boat can go back on the same path in the same amount of time.

Furthermore, when V(x) is symmetric wrt x=1/2 i.e. V(x)=V(1-x) (as in the provided parabolic example) then it seems that  y(x) = - y(1-x).

[In the worksheet, from sin(psi-alpha)=k ==> psi-alpha = arcsin(k) OR  psi-alpha = Pi - arcsin(k) ...].

[Edited]

@Rouben Rostamian  

Excellent presentation. Congratulations!

A bare definition appears here.