Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

Posting a picture of your equations is not very productive.  I have a couple of ideas toward the solution but I am surely not going to retype your equations.

Consider uploading a maple worksheet so that people can work with it.  Look for the large green up-arrow in the screen where you post/edit your question.

In the mean time, you should also think about the following:

  1.  You say the derivatives of u are small.  What about u itself?  Is it small as well?  (I am writing u for your uu0 for simplicity.)
  2.  You say beta is small.  Are the derivatives of beta small as well?  Note that a function may be small but its derivative may be quite large;  think of y(x) = 1/k *  sin(k^2*x) when k is large.
  3. How does the smallness of beta compare to the smallness of u?  Suppose the smallness is measured by a small number, let's say h.  Compare the case when
        u = u0 + h*u1,   beta = beta0 + h*beta1,
    versus
        u = u0 + h*u1,   beta = beta0 + h^2*beta1.
    These will lead to two different linearizations in general.

@acer I did not know about Explore, I am sure it will come in handy for me in the future.  That's the sort of thing I had in mind when I referred to Maplets.

@Kitonum That's nice.  I should remember to say "uses" instead of with(...) in my procs.

@brian bovril Suppose that a point moves on the plane so that its polar coordinates are (r(t)),θ(t)).  Then its Cartesian coordinates are x = r cos θ, y = r sin θ, whence its velocity components are

  x' = r' cos θ - r θ' sin θ,
  y' = r' sin θ + r θ' cos θ.

were x' means dx/dt, etc.  Squaring and adding up we get x'^2 + y'^2 = r'^2 + r^2 θ'^2.

As to our pirate ship, we have x'^2 + y'^2 = V^2, and  r(t) = v t. therefore  V^2 = v^2 + v^2  t^2 θ'^2, as asserted.

Regarding animation, Kitonum and acer have posted nice responses although those do not quite parallel the mathematical details of the solution that I had posted.   Here is how I did my animation.  (I have added markers for the ships following Kitonum's idea.)  This one may be easier to follow if you want to match with my explanation.

restart;
with(plots):
pirate_pursuit := proc(V, v, L, a, nframes)
  local b, T, t_capture, P_r, P_a, M_r, M_a, frames;

  T := L/(V+v);
  b := t -> sqrt((V/v)^2-1) * ln((V+v)/L*t);
  t_capture := solve(b(t)=a, t);

  # pirate's polar radius and polar angle as functions of time
  P_r := unapply(piecewise(t < T, L - V*t, v*t), t);
  P_a := unapply(piecewise(t < T, 0, b(t)), t);

  # merchant's polar radius and polar angle as functions of time
  M_r := t -> v*t;
  M_a := t -> a;

  # animation frames
  frames := seq(
    display([
      plot([P_r(s), P_a(s), s=0..t], coords=polar, color=red, thickness=5),
      plot([M_r(s), M_a(s), s=0..t], coords=polar, color=blue, thickness=5),
      pointplot([[P_r(t), P_a(t)], [M_r(t), M_a(t)]],
        coords=polar, color=[red,blue ], symbol=solidcircle, symbolsize=30)
    ]) , t=0.. evalf(t_capture), evalf(t_capture/(nframes-1)));

  # animate the paths of the two ships:
  display([frames], insequence, scaling=constrained, axes=normal);

end proc:

# do a demo:
pirate_pursuit(3, 1, 1, 2*Pi-Pi/4, 50);

@gkokovidis The article you have cited is unrelated to the problem at hand.  That article assumes that the persuer can see the persued at all times.  In our problem the persued disappears in the fog.  A different strategy is required.

@wolverine Yes, that is a terrible mistake -- there are no t1 and t2 in your worksheet.

Since you haven't used any words to explain what it is that you want, I am attempting a guess here in total darknes.

restart;

de := diff(y(t),t,t) - e*(1-y(t)^2)*diff(y(t),t) + y(t) = 0;

de1 := subs(y(t) = eta[0](t) + e*eta[1](t), de);

de2 := collect(de1, e);

coeff(lhs(de2), e, 0) = 0;

coeff(lhs(de2), e, 1) = 0;

@Markiyan Hirnyk Acer's post shows much more than just Maple skills.  There is a darn good deal of mathematical inventiveness there.  Same goes for Axel's approach.

@patient I am a bit suspicious of the solution in your worksheet.  There you set b[] equal to the initial condition, and b[] does not change during the rest of the computation.  Do you really mean that?  That looks very odd to me.

Aside: Your calculations may be simplified quite a bit if you compute the matrix P explicitly, as I outlined earlier.  Since P is constant, its inverse, let's call it Pinv, can be calculated just once, as in Pinv := P^(-1), before the iteration begins.  Then the main loop of the computation takes the form

    set U = a  # Initial condition   
    p := plot(<X|U>;
    for i from 1 to N do
      Q := RightHandSide(U);   # RightHandSide() produces the right hand side of the equation
      U := Pinv.Q;
      p := p, plot(<X|U>;
    end do;

Upon exit from that loop, p will contain a sequence of plots.  To see an animation, do

    plots:-display([p], insequence=true);

@Preben Alsholm That's very clever.  I believe that this is just the right way of solving the problem although I don't quite understand the solution's method.  In order to gain a better understanding, I attempted to modify it slighltly in a way that is more comprehensible to me, but I got answers which were obviously wrong.  I will think about this some more and ask you for help if needed.

@Carl Love  At first reading I thought that your comment was tongue-in-cheek, but then I realized that perhaps you are bring serious.  Like you, I haven't had a need for programmatic use of the with command.  I do, however, find with's output in case of a less-frequently used package quite handy.  I know, for example, that the Laplace transform is in the inttrans package. Now, is it "laplace" or "Laplace"?  And what is the command for the inverse Laplace transform?

Another use of with's command is for classroom teaching.  I don't expect an undergraduate student who is taking four other courses besides mine to remember why I load the plots package when doing a demo for them.  Letting with(plots)'s output to spill onto the screen gives me the opportunity to remind the class that this is so that we may access the display and pointplot commands because they may not remember that from the previous demo from two weeks ago.

 

@Preben Alsholm I will address each of the issues that you have brought up.

1. By de2, v cannot be zero on any interval no matter what value is assigned to signum(0).

You are right.  That's why I said it remains a mystery to me how Maple gets that solution.

But as I also explained, Maple's solution agrees with the hand-calculated exact soluition  of the problem (not the solution of the DE) which I found by splitting the motion into a sequence of stick and slip phases.  The exact solution for the absolute displacement of the house, that is x(t)+X(t), is a piecewise function consisting of sinosoidal and parabolic patches.

2. If you look at an interval where v appears to be zero ... we see rapid oscillation.

Those oscillations are numerical artifacts.  Their amplitude is roughly 0.004, which is in effect zero considering the choice of stepsize=1e-3.  If we set stepsize=1e-4, the amplitude is cut down by a factor of 10.  And if we set stepsize=1e-5, the amplitude is cut down by a factor of 100.

We see that Maple is attempting to produce an identically zero function v(t) on the interval 3.25..3.9.  As you have observed, an identically zero function is inconsistent with the differential equation de2.  But mysteriously, it is consistent with the exact solution which one obtains by hand by splitting the motion into a sequence of stick and slip phases.

3. Making changes in order to match the original problem I tried using g = 9.81, mu=0.8, and cos instead of sin...The results where disturbing.

The results are just right.  With your choice of parameters there will be no slippage, and therefore the solution v(t) should be identically zero.   Maple shows v(t) oscillating with an amplitude of roughly 0.0002, which is in effect zero plus numerical noise.

To elaborate, let f be the horizontal force exerted by the ground on the house.  Newton's law says f = m*(a + A), where as before, A is the acceleration of the ground, and a is the house's acceleration relative to the ground.  In particular, if the house does not slide, then f = m*A.  The house will begin to slide only when |m*A| exceeds the force of friction, which is  mu*m*g.  With your choice of parameters, |m*A(t)| is always less than mu*m*g, therefore the house never slides.

4. Thus I don't trust the results. And I should add: The mathematical model consisting of just de2 (and de1).

I agree.  de2 is definitely inadequate becasue as you have observed, it is inconsistent with v(t) being zero on an interval.  I do not understand how Maple "solves" that incorrect model and magically obtains the  correct solution  solution to the correct model.

 

@Preben Alsholm Your criticism is justified.  The problem is more complex than I presented it.  For some magical reason, however, the following naive approach produces what appears to be the correct solution.  I will explain further down why I believe the solution is correct.

restart;
with(plots):  with(plottools):
de1 := diff(x(t),t) = v(t);
de2 := m*diff(v(t),t) = -mu*m*g*signum(v(t)) - m*A(t);
X := t -> P*sin(omega*t);
V := D(X);
A := D(V);
ic := x(0)=0, v(0)=0;
m:=1: g:=1: P:=1: mu:=0.7: omega:=1:
dsol := dsolve({de1,de2,ic}, numeric, stepsize=1e-3, maxfun=0);
odeplot(dsol, [t, x(t)], t=0..6*Pi/omega);

odeplot(dsol, [t, v(t)], t=0..6*Pi/omega);

 

These pictures agree with the x(t) and v(t) graphs given earlier in this discussion.

Furthermore, I also obtained these same graphs through a totally different approach which I think is quite defensible.  I split the motiun into distinct phases, distinguishing between the case where the house sticks to the ground, and the case when it slips against it.  In the former case the house's moiton coincides with that of the gound,  In the latter case the house moves under the action of the contact force -mu*m*g*signum(v), which is a constant, therefore the velocity is linear in time, and the displacement is quadratic in time.   Therefore the overall displacement over time is a piecewise curve conising of sinusoidal and parabolic patches.  The resulting graphs coincide with those shown above.  It remains a mystery to me why the Maple script shown above produces the same result.

Here is what the motion looks like:

@Preben Alsholm I analyzed the problem from scratch and found the correct equations of motion which I will now describe.  I don't know how to solve the resulting equations in Maple, but probably you know a trick.

Referring to the diagram posted earlier, let X(t) be the displacement of the ground, and x(t) be the displacement of the house relative to the ground.  [Thus, I am replacing the diagram's x_sol(t) with X(t).]  We have
 
X := t -> P*cos(omega*t);   # P for amplitude because we have a better use for A

Let V := dX/dt, A := dV/dt, v := dx/dt, a := dv/dt.

The equations of motion are:

de1 := diff(x(t),t) = v(t);
de2 := diff(v(t),t) =
        -mu*g*sigum(v(t)) - A(t),  if abs(A+a) < mu*g
        0                          otherwise

The first part of de2 is almost like that shown in this discussion's initial post, except that the rightmost term seems to have the wrong sign there.

@Oloyemike I understand your model a little better now, but things are still unclear.  At places you have Q(r, x, tau), at other places you have Q(r, tau), and yet at some other places you have just Q.  Shouldn't there always be three arguments to each Q?

Similarly, at places you have U(x,tau),  and at some other places you have just U.  Shouldn't there always be two arguments to each U?

Mathematics is a precision instrument.  Sloppiness is not acceptable.

Furthermore, based on what you have shown, I can tell you the solution of your system of PDEs, and I don't need Maple to do it.  The solution is U(x,tau) = 0 and Q(r,x,tau) = 0.  Just substitute these in the equations to verify. If this is not what you had expected, then your model needs further adjustments.  Oops, sorry, I take that back -- U(x,tau) = 0 and Q(r,x,tau) = 0 is not a solution.

 

 

 

@Carl Love Let's say we have a ball of radius R, and let r be the distance from the center of the ball.  As we move from the center to the ball's boundary, r goes from 0 to R.  Let's introduce rho = r/R.  Then as we move from the center to the boundary, rho goes from 0 to 1.  In this context, rho is called the dimensionless radius.  Similarly one defines dimensionless time, dimensionless force, etc.  The terminology is quite common in engineering and physics.

First 89 90 91 92 93 94 95 Last Page 91 of 99