@oceanxinlie You asked:

Can I know where do you find the elementary function for solving this differential equation, is it very famous?

Your differential equation y'' = a*y^q has quite a simple form, therefore is not difficult to guess that it must have a solution of the form y(x) = A*x^B.

To determine A and B, we substitute that guess into the equation and equate the exponents of x on both sides and arrive at B = 2/(1-q).  Then we equate the coefficients of x on both sides, whence we arrive at the value of A that you see in the worksheet.

The solution y(x)=A*x^B has the property that y(0)=0 and y'(0)=0.  This does not solve your boundary value problem because you want y(0)=0 and y'(L)=b.  But that's alright.  You want L=2945 but let's ignore that for the moment.  Let's find a different value of L, let's call it Lcrit, so that y'(Lcrit)=b.  Since y(x)=A*x^B, this says A*B*Lcrit^(B-1) = b.  We solve this equation for Lcrit and arrive at the value of Lcrit given in the worksheet.

This would be the end of our calculations if your L were exactly equal to Lcrit.  But it isn't.  Your L is quite a bit larger than Lcrit.  But that's not a problem.  It is not hard to see that if a function y(x) solves the differential equation, then the shifted function y(x-k) also solves the equation for any k.  So we pick k = L − Lcrit so that the shifted function satisfies your desired boundary condition when x = L.

We fill the gap on the x axis between the origin and the shifted function by zero.  This is a well-known idea which Preben immediately recognized with his comment about y(x0) and y'(x0).

You also asked:

Can you show me the code of numerical solution?

That's exactly what Preben did in his first example.  He took L=450 which is less than Lcrit, and found the solution numerically.

The usual numerical methods don't work when L is larger that Lcrit, but that's not a problem because in that case we have an exact solution as shown in the worksheet.

In the rest of your comments you say that you had hoped for a solution which is not flat zero on a part of the interval.  As far as I can see, if L > Lcrit, that zero part cannot be avoided.  Mathematically there is nothing wrong with that.  If you are applying the result to some physical model, you will have to see what that is saying about the model.

@mohammad1369 If I understand your question correctly, you want to write A_x instead of A.  You enter A_x in Maple by typing A__x, that is, A followed by two underscore characters and an x.  For example:

A__x := <a, b; c, d>;

@Les Sure, there is nothing preventing from putting that construction inside a procedure.

@dpaddy Before you get busy coding your function in C, see if this can help you:

h := t -> exp(I*t);
plot(Re(h(t)), t=0..2*Pi);

@Earl If the disk's initial angular velocity is small, then it will rock back-and-forth, rather than take off and roll forward forever.  Note that unlike your drawing, I am taking the blue disk to be at the lowest position at time zero.

The transition from rocking to rolling motion is obtained by giving a sufficiently large initial angular velocity that enables the blue disk to rise from the bottom position to the top position.  At the borderline case, the blue disk will creep to the top position and stop there.  (By the way, in your post where you write Pi/2 it should be Pi.)  That borderline case is determined by equating the initial kinetic energy to the change in potential energy.  The expressions for the two energies are available in the worksheet.  Equating the two and solving for omega^2 we get

omega^2 = 8*g*m*d/((2*R^2-4*R*d+2*d^2+r^2)*m+3*M*R^2), that is:

Plugging in the parameters used in the worksheet we get omega^2 = 0.4232804233, and therefore omega = 0.6506000486, which is what you have observed.

As to g=1, that depends on the choice of units.  The value of g may be 9.8, or 980, or 32.2,, or anything else, depending on the choices of the units of mass, length, and time.  Since we are not saying what those units are, there is no reason to be so particular about the value of g.

Regarding your and JohnS's question about what happens after the vertical contact force hits zero, I expect the disk to become airborne and spin freely for a while, and then fall back and regain contact.  I haven't done the calculations for that.  I may look at it this weekend if I find the time.

@max125 The command int~ is  an instance of Maple's element-wise operation.  Look up the help page on "Element-wise Operators".  In short, int~ applies the integration operator to the individual components of a vector.

I must add that vv's answer is the best solution to your question.  That, too, uses int~ to integrate the vector.

Your solution looks good to me.  Here's a slight variation:

restart;
a := < 4*t, sin(t), cos(2*t) >;
v0 := <1,0,0>;  r0 := <0,1,0>;
# calculate v:
int~(a, t) + C;
eval(%, t=0) = v0;
simplify(isolate(%, C));
v := eval(%%%, %);
# calculate r:
int~(v, t) + C;
eval(%, t=0) = r0;
simplify(isolate(%, C));
r := eval(%%%, %);

There is nothing wrong, mathematically, with your equations.  I assume you are entering them incorrectly.  Post your worksheet so that it can be examined.

Click on the big fat green arrow in the window where you compose your message in order to upload the worksheet.

@vv Consider the parameters m=k=g=mu=1.  Let x(0)=1/2, x'(0)=0.  What is the solution?

From the physics of the problem we can tell that the solution should be x(t)=1/2 for all t > 0.  But in what sense does this solution satisfy the differential equation if signum(0) is undefined?

As has been observed by Preben and others, the inequality containing the second derivative of x(t) in the original post is incorrect.  It is supposed to be the first derivative.

This leads to the differential equation

m*diff(x(t),t,t) = - k*x(t) - mu*m*g*signum(v(t));

where v(t) = diff(x(t),t).

That differential equation does not reflect the problem's physics due to the arbitrariness in the definition of the signum function when v=0.  That's a very significant omission.

The proper formulation of the mathematical model calls for an extension of the idea of differential equations to objects called differential inclusions.  Here is a quick summary:

Define a set-valued extension of the signum function, let's call is Signum (with capital S) through:

Signum(v) = {+1} if v > 0,
          = {-1} if v < 0,
          = [-1,1] if v = 0.       #  [-1,1] is the closed interval from -1 to +1.

Now that Signum is a set-valued, the differential equation changes to a differential inclusion:

m*diff(x(t),t,t) ∈ - k*x(t) - mu*m*g*Signum(v(t));

There exists extensive literature on the theory and numerical algorithms for solutions of differential inclusion.  A good place to start reading about this is the article

D. E. Stewart and J. C. Trinkle, Dynamics, friction, and complementarity problems

published in Complementarity and variational problems, SIAM, 1997.

A PDF of the mentioned article can be found on the Internet for download.

@Earl To gain some intuitive insight into the process, let's consider the case where the water's velocity is constant throughout, that is, it does not vary from point to point.  If the boat attempts to go to the point directly across the river, there is no gain in heading downstream at any time because then it will have to exert extra effort to bring itself back.  So we expect the heading to be toward upstream throughout the travel.

Let V be the (constant) vector of the water's speed and S the vector of the boat's speed relative to the water.  The vector V points along the y axis.  The vector S has a component along the y axis to counteract V, and a component along the x axis to carry the boat across.  The y components of S and V have to cancel each other to prevent the boat from drifting downstream, therefore the length of S should be longer than the length of V.  That is, the boat's speed should be greater than the water's speed.

That is the intuitive explanation of your observation that one cannot reduce the boat's speed too much.

The requirement of sufficiently large S becomes more severe if the destination is upstream because the boat needs to compensate for the water's speed, and then some, in order to propel itself upstream.

@vv I don't understand your comment about non-uniqueness and the statment that "because the boat can go back on the same path in the same amount of time".  If the boat goes back on the same path, the the path is still y(x), not y(1-x).

Your observation about the symmetries of V(x) and y(x) is absolutely true.  This may be verified directly in Maple with the help of PDEtoos:-dchange.  To simplify the calculations, and without a loss of generality, it helps to shift the x axis so that the y axis is the axis of symmetry, and assume that V(x)=V(-x).  Then we introduce  new variables x* = -x and y*(x*) = -y(x), and rewrite the differential equation in terms of the new variables. Then we verify that the differential equation is invariant under that transformation.  The attached worksheet, which is a continuation of the previous one, has the details.

boat2-ver3.mw

Earl and vv (and anyone else reading this thread):

There is an interesting consequence to the analysis of the optimal path.  Suppose the water's velocity is uniform throughout the river, that is, the water's velocity vector is <0,V> with a constant V everywhere.  Then it turns out that the fastest path to go from a point on one bank to any point on the opposite bank (not necessarily directly across) is a straight line.

In particular, the fastest path to the point directly across the river is the straight line which is perpendicular to the stream.  Initially I found that to be counter-intuitive but after some thinking about it, I see how it works.

This assumes that the boat's speed, S, is faster than the water's speed, V, otherwise there is no way to sail to the point directly across.

A couple of things:

1. The worksheet is not accessible.  See if you can delete the link and upload it again.
2. In the problem's statement we have "...crosses the river heading at all times towards its destination".  We also have "captain could constantly vary his heading so as...".  It seems to me that these two statements regarding heading are inconsistent.

The question that you have posted is incomplete:

1. It's hard to tell what the truncation error is if you don't say what it is that is being truncated.
2. The meaning of the expression that you have shown is unclear.  What does it represent? Shouldn't there be an equal sign somehwere?
3. I can guess what h is, but what are f, g, and y?
4. Where do all those fractional coefficients come from?