Preben Alsholm

13728 Reputation

22 Badges

20 years, 253 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@krw2015 I tried the lines

with(DifferentialGeometry); with(Tensor); with(Tools);

DGsetup([t, r, theta, phi], M1, verbose);

in Maple 17.02 and in Maple 2015.1. I didn't receive any errors in either version.

@ It could be pointed out that if we don't restrict ourselves to the interval -Pi/2..Pi/2 but consider solutions valid for all t then the problem
diff(y(t),t)=z(t), y(t)^2+z(t)^2=1 with y(0)=0,z(0)=1
has infinitely many solutions with y continuously differentiable and z continuous.
An example:

restart;
ode:=diff(y(t),t)=z(t);
eq:=y(t)^2+z(t)^2=1;
ics:={y(0)=0,z(0)=1};
dsolve({ode,eq});
Y:=piecewise(t<Pi/2,sin(t),t<4,1,t<4+Pi,sin(t+Pi/2-4),-1);
Z:=diff(Y,t);
plot(Y,t=0..10);
plot(Z,t=0..10,thickness=3);

#########Added note ########
It is interesting to see what dsolve does to the problem if we solve eq for z(t) and choose the nonnegative solution:
ode2:=diff(y(t),t)=sqrt(1-y(t)^2);
dsolve({ode2,y(0)=0});
#We get y(t) = sin(t), which is correct on the interval -Pi/2..Pi/2, but not on any larger interval containing t=0.
#Now numeric solution:
res:=dsolve({ode2,y(0)=0},numeric);
plots:-odeplot(res,[t,y(t)],-5..5); #The solution is constant outside the interval -Pi/2..Pi/2 as one would expect.




@Arthurok What matters is that f1(y) depends on y only and that f2(x) depends on x only.
The equations eq1 and eq2 gives you the only possible solutions for u(x,y) and v(x,y).

Inserting those into pde3 and differentiating the resulting equation w.r.t x and y as I did above you get the equation
-(2/3)*(-72*nu*y+336*y)/h^3+48*y*nu/h^3 = (48*(1+nu))*y/h^3
which must be satisfied for (at least) a range of y-values. Notice that f1(y) and f2(x) are gone.

Now the difference between the right and left hand sides of this equation should then be identically zero (on an interval at least). But the difference is (as found above)
-16*y*(-17+3*nu)/h^3

This is identically zero if and only if nu = 17/3. So in general the problem has no solution.

If on the other hand nu=17/3 then pdsolve can find the solutions:
res:=pdsolve(eval({pde1,pde2,pde3},nu=17/3));

However, your initial value problem has no solution:
pdsolve(eval({pde1,pde2,pde3,u(L,0) = 0,v(L,0) = 0,D[1](u)(L,0) = 0},nu=17/3));

Notice that the * you had in the third condition makes no sense. The derivative of u w.r.t. the first variable (x) at (x,y) = (L,0) is D[1](u)(L,0).





@Rouben Rostamian  I would like to call attention to the fact that I have edited the code for the two events. The changes are described in the edited version. One is logical the other is numerical.

I just tried copy and paste on your list of lists and replaced the comma at the very end with a semicolon.
I got one execution gropup, but as you say several prompts. Those latter prompts don't matter. It executed fine in Maple2015, worksheet mode and 1d input.

@Rouben Rostamian  First of all thank you very much for clearing this up.
I think I have a solution now using 'events'.
For convenience in changing parameters (e.g. g or P) I also use the parameters option, but that is for convenience only and can be left out. If so parameters must be assigned ahead of using dsolve.
#Edited relative to the original version:
1. Changed the second event from [a/mu/g=-1..1,[b(t)=1,n(t)=signum(v(t))]] to [a/mu/g=-1..1,[b(t)=1,n(t)=signum(a)]], because when abs(a)>mu*g the velocity will increase from its zero value if a>0 and decrease from zero if a < 0. Thus n(t) will have the value signum(v(t)) right after the event. Logically this makes sense. The success of the original could be called accidental.
2. Changed the first event from [v(t)=0,[b(t)=0,n(t)=0]] to [v(t)=0,[b(t)=0,n(t)=0,v(t)=0]] to ensure a perfect zero for v(t). This is done for numerical reasons.


restart;
X := t -> P*cos(omega*t);   # P for amplitude
#Let V := dX/dt, A := dV/dt, v := dx/dt, a := dv/dt.
#The equations of motion are de1 and de2:
de1 := diff(x(t),t) = v(t);
A:=diff(X(t),t,t);
a:=-mu*g*n(t)-A; #For n(t) see the dsolve command
de2:=diff(v(t),t)=b(t)*a; #For b(t) see the dsolve command
##
res:=dsolve(
       {de1,de2,x(0)=0.5,v(0)=0.25,b(0)=1,n(0)=1},numeric,
       discrete_variables=[b(t)::boolean,n(t)::integer],
       events=[[v(t)=0,[b(t)=0,n(t)=0,v(t)=0]],[a/mu/g=-1..1,[b(t)=1,n(t)=signum(a)]]],
       parameters=[P,omega,mu,g]
       );
res(parameters=[P=1,omega=1,mu=0.8,g=1]); #setting parameters
plots:-odeplot(res,[t,v(t)],0..15,thickness=3);
plots:-odeplot(res,[t,x(t)],0..15,thickness=3);
plots:-display(Array([plots:-odeplot(res,[t,n(t)],0..15,thickness=3),plots:-odeplot(res,[t,b(t)],0..15,thickness=3)]));
##The graph for v(t) with the parameters set as above:



###########################
#With parameters = [P=1,omega=1,mu=0.8,g=9.81] we get what you described: The house comes to relative rest and remains at relative rest.
#Now changing also P:
res(parameters=[P=10,omega=1,mu=0.8,g=9.81]);
plots:-odeplot(res,[t,v(t)],0..15,thickness=3);
plots:-odeplot(res,[t,x(t)],0..15,thickness=3);
plots:-display(Array([plots:-odeplot(res,[t,n(t)],0..15,thickness=3),plots:-odeplot(res,[t,b(t)],0..15,thickness=3)]));


@Rouben Rostamian  By de2 v cannot be zero on any interval no matter what value is assigned to signum(0).
If you look at an interval where v appears to be zero as it does on 3.25..3.9:
odeplot(dsol, [t, v(t)], t=3.25..3.9);

then we see rapid oscillation. Also the following plot shows that.
odeplot(dsol,[t,signum(v(t))],t=3.25..3.9,thickness=3,style=point);

But as I understand the physics the house must be at relative rest (i.e. v(t) = 0) on intervals.
#########################
Making changes in order to match the original problem I tried using g = 9.81, mu=0.8, and cos instead of sin:

X := t -> P*cos(omega*t);
V := D(X);
A := D(V);
ic := x(0)=0, v(0)=0;
m:=1: g:=9.81: P:=1: mu:=0.8: omega:=1:

The results where disturbing. This one in particular:
odeplot(dsol, [t, v(t)], t=0..6*Pi/omega);

Another observation: While the results when using stepsize and abserr looked much like each other for your choice of parameters, they didn't for the choice above. The graphs of x(t) were rather different, while both graphs of v(t) were rather wild (but different).

Thus I don't trust the results.

And I should add: The mathematical model consisting of just de2 (and de1).






@phil76600 I tested Rouben's code in Maple 12.02. It didn't complain about 'stepsize'. Which Maple do you have 12.01 or 12.0?

But try replacing the word 'stepsize' by 'abserr'.

maxfun = 0 means that there is no limit to the number of function evaluations when computing the solution. You could say that a more intuitive way to do this would have been maxfun = infinity.

@Rouben Rostamian  This is very interesting. I shall have a closer look.

I was surprised by the effect of giving the stepsize. I didn't know that option existed except for the classical methods!

@phil76600 Your revised problem:

ode1 := m*diff(x(t), t, t) = k*Xr-mu*m*g if mu*m*g<k*Xr else diff(x(t),t) = V(t) = 0

(if I understand it correctly) can be solved by using 'events'.
See the help for that.
?dsolve,events

The event [k*Xr=mu*m*g..infinity,diff(x(t),t)=0] means that if kXr falls outside of the interval mu*m*g..infinity at some time t1 then the action is to set diff(x(t),t) = 0 at that time t1. 

restart;
Xr := U*t-x(t);
k := 10; m := 1; g := 10; mu := .2; U:=1;
ode1 := m*diff(x(t), t, t) = k*Xr-mu*m*g;
res:=dsolve({ode1,x(0)=0,D(x)(0)=0},numeric,events=[[k*Xr=mu*m*g..infinity,diff(x(t),t)=0]]);
plots:-odeplot(res,[t,x(t)],0..15);
plots:-odeplot(res,[t,diff(x(t),t)],0..15);
plots:-odeplot(res,[t,mu*m*g-Xr],0..15,thickness=3);




@phil76600 Let us not mix the two discussions, but let me just point out that mu*m*g>m*ω^2*A only depends on the given parameters, so either this relation holds under the whole process or it doesn't.

@phil76600 Actually they are not the same. Take a look at the of the interval t=0..15 (e.g.).
I saved both plots and displayed them together:

Notice that for mu = 1000 the piecewise gives 0 for all t in the interval 0..15. This is not the case for mu=0.2.
Using
pw:=piecewise(mu*m*g<k*Xr,mu*m*g,0);
you can try
plots:-odeplot(res,[t,pw],0..15,thickness=3);
You will notice that for mu = 1000 the values are all zero (therefore I added thickness=3). For mu=0.2 this is definitely not the case.

@Rouben Rostamian  Using your notation a = diff(v(t),t), so from your equation de2 we get a+A = -mu*g*signum(v(t)). Thus abs(A+a) = mu*g*abs(signum(t)), which is equal to mu*g if v(t)<>0 and equal to 0 if v(t)=0 if we accept that signum(0) = 0.
Thus the first part of the piecewise right hand side only holds if v(t) = 0, whereas the second case (the otherwise,output 0) holds whenever v(t)<>0.
Surely, you didn't mean it like that.

@tobbie247 If you my code from above and try
res(0);

you will get very nearly the values for A and B you gave in a reply to Carl:
0.742408087440 and 0.943866213084

I don't know anything about the Differential Transform Method either.

@phil76600
I'm not happy with my solution either. I'll have to rethink this.

## Why do you have cos(t)^2 in the Maple code when cos(omega*t) is not squared in the image of the formula?
## Are you sure mu = 0.8 is correct?
## Is the model actually correct as stated (apart from the square mentioned above?)
As I understand it presently, the ground on which the house is standing is moving in a periodic way as in
xsol(t) = A*cos(omega*t).
x(t) is the position of the house relative to the ground.
Thus the abolute position of the house is X(t) = x(t) + xsol(t).
Using Newton's 2nd law we get
m*X'' = friction
so that
m*(x'' + xsol'') = friction
Thus
x'' = friction/m - xsol'' = friction/m + A*omega^2*cos(omega*t) .
The problem for me is what to set friction to be when there is no relative motion, i.e. when v(t) = 0.
Am I misunderstanding something?

First 112 113 114 115 116 117 118 Last Page 114 of 230