Preben Alsholm

13728 Reputation

22 Badges

20 years, 254 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Christopher2222 Sorry about not catching the obvious mistake right away. Assuming that up is the same as up on your figure in the worksheet, then the altitude is given by h=R*sin(theta(t)) when theta is measured from the positive x-axis as you indicate in your figure.

The potential energy is then m*g*h + C, where the constant C can be chosen arbitrarily.

Had the angle theta been measured from the positive y-axis then the altitude would be given by h=R*cos(theta(t)) and the potential energy would again be m*g*h + C. Keeping the initial conditions then the events trigger would be theta(t)=Pi/2 (and theta(t)=-Pi/2, but that won't ever happen for your initial condition).

@patient Clearly the initial conditions for S and U must be different from the initial conditions for s and u. The change of variable from z to x obviously has an effect.
I won't change your (to my mind) unfortunate name x0 for the initial z-value (z0 appears to me to be more natural).
I shall call the x-value corresponding to z=x0 for x00.
That value is found from x^2=t*z to be x00 = sqrt(t*x0) if we choose x >0, which I shall do.

The derivative of S at x = x00 is found by using the chain rule to be
S'(x) = s'(z)*2*x/t
thus
S'(x00)= s'(x0)*2*x00/t.

Now we notice that the parameter t appears in the initial conditions as S(sqrt(t*x0)) = ...
It turns out that dsolve doesn't accept that, whereas the form S(x00) = ... with x00 as a parameter is accepted. Whether that should be considered a bug or not is irrelevant right now: we can work around it by using x00 as parameter instead of t.

Thus we can do:
SYS2:=subs(t=x00^2/x0,SYS);
ICS:={S(x00)=exp(x0^2-16),U(x00)=9*exp(x0^2-10),D(S)(x00)=-4*exp(-12)*2*x0/x00};
res:=dsolve(SYS2 union ICS, numeric, method = rkf45_dae, parameters = [x00]);
res(parameters=[sqrt(x0*(-.1))]); # This means that t = -0.1
plots:-odeplot(res,[[x,100*S(x)],[x,U(x)]],0..10);

Notice that I have multiplied S by 100 so it doesn't get drowned by U.

The procedure p can be modified as follows:
p := proc (t)
   res(parameters = [sqrt(t*x0)]);
   plots:-odeplot(res, [x, U(x)], 0 .. 4)
end proc;

plots:-animate(p, [t], t = -.1 .. -0.1e-1,view=0..0.05);


@Christopher2222 The zero for potential energy can be set arbitrarily. It is of no consequence for the equations of motion. In fact just look at the case in hand. L = T-U. You only use two derivatives of L. Both would eliminate an additive constant to U. The main thing is that since the gravitational pull is downwards, U should increase with the altitude.

@Christopher2222 Obviously the potential energy should increase with the altitude. Without a change of sign it doesn't. Whether U is positive or negative is irrelevant.

As Rouben points out, there is nothing wrong.

I can produce the same error by moving the local declaration:

game := proc()
  #local player1, player2, roll;
  roll := rand( 1..6 );
  local player1, player2, roll;
  player1 := roll():
  player2 := roll(2):
  if player1>player2 then "A wins"
  elif player1=player2 then "Tie"
  else "B wins"
  end if;
end proc:

I think you ought to reveal some more of what you are doing. I sure don't understand it.

@patient Changing variable can be done with dchange from the PDEtools package:

SYS:=PDEtools:-dchange({z=x^2*1/t,u(z)=U(x),v(z)=V(x)},sys,[x,U,V]);
#I just check that the DAE method will work:
solve({SYS[1],diff(SYS[2],x)},{diff(V(x),x,x),diff(U(x),x)});
#Remember that the initial conditions have to satisfy this equation:
SYS[2];
ICS:={U(0.1)=1,V(0.1)=1,D(V)(0.1)=0}; #Just an example.
res:=dsolve(SYS union ics, numeric,method=rkf45_dae,parameters=[t]);
res(parameters=[-1]); #Setting paramer t to -1
plots:-odeplot(res,[[x,U(x)],[x,V(x)]],0.001..3);
#Experimenting with animation in t with the given initial conditions:
p:=proc(t) res(parameters=[t]); plots:-odeplot(res,[[x,U(x)],[x,V(x)]],0.001..3) end proc;
plots:-animate(p,[t],t=-1..-0.01);

@emmantop OK. The reason I asked was that D has a very special meaning in Maple (as you know: you used it in your worksheet). So it is better to use some other symbol for the sake of clarity.

@escorpsy There is no difference between diff(A(r),r)^2 and (diff(A(r),r))^2 . In other words, the outer parenthesis in the latter is superfluous.

Your system PDE only has derivatives w.r.t. y, thus it is an ode system.

Your IBCs contains the equation

 (D[1](phi))(R(z)/R[0]) = R[0]*delta/(D*phi[w])

What is (D*phi[w]) meant to be?

@escorpsy Maybe I misunderstood the meaning of your term "test solution". I thought that you knew (or thought) that the functions given in puta made up a solution to your ode.
If puta was a solution then odetest should return {0}.

I fail to see the connection between the original ode and the rest of the code.

Your first order system would be
y'(t) = u(t)
u'(t) = -3*u(t)-2*y(t) + cos(t)

Just one indication of the missing connection: Where is cos(t) in the code?

@surnizai The range was chosen to include the result of fsolve(x*y1,x) which I actually did before I plotted.

From the plot I could see that there was (at least) one more zero and that one was near x=1+I*epsilon. So giving x=1 as a start to fsolve might get us the second roots. It did!

@patient I noticed that you changed both of the equations so you have to think again. I shall be busy with travelling the next few days so won't have time to look at it.

@matmxhu I tried the following. When doing it several times the error was successfully caught some of the time and sometimes not:

restart;
g:=proc() solve(exp(x)=1,x) end proc:
try
 timelimit(0.02, g());
catch "time expired":
 123456789
end try;

I don't know why the error is not caught all the time.
This was done in Maple 2015.
I tried in Maple 15 several times. The error was caught every time I tried.
In Maple 16, 17, and 18 the error was caught sometimes and sometimes not.

With the following code the error was caught all the times I tried in Maple 15, 16, 17, 18, and in 2015.

restart;
f:=proc() local i; for i to 10000 do evalf(exp(i)) end do end proc:
try
 timelimit(0.02, f());
catch "time expired":
 123456789
end try;

First 120 121 122 123 124 125 126 Last Page 122 of 230