Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@ewas OK, but you still need to find the four extra boundary conditions relevant to your problem.

There are som missing parameter values. In particular B1 and C1 since they are rather crucial being the coeffcients of highest order. Thus it makes an enormous difference if they are zero or not. On the other hand the information that C2=0 is irrelevant since it doesn't appear in your system.
The spelling is 'dsolve', not 'Dsolve'.
What are the values of az, aw (you wrote 'and'), and ag?
Two of the boundary conditions could be
az * Tz (0) + aw * Tw (0) + ag * Tg (0)=-10, az * Tz (20) + aw * Tw (20) + ag * Tg (20)=10.
You need 4 additional ones.

I don't understand the do-loop.

Since tan(x)  appears in your system and has a singularity at all odd multiples of Pi/2 you should reconsider using the interval x=0..20.

@maple fan I had no problem plotting the solution on the range 0..8.
The plot was rather boring: Very soon the graph became linear with slope 0.1.
What did you expect should happen?
What has mu got to do with the problem?

The following code does the same thing:
restart;
u:=1003-1000*x(t)-30*diff(x(t), t)-25*signum(diff(x(t), t)-.1)-.3*signum(diff(x(t), t))*exp(-2*abs(diff(x(t), t)));
sys := diff(x(t), t, t) = (1-b(t))*u, x(0)=1, D(x)(0) = 0;
stick := [[diff(x(t), t) = .1, abs(1-x(t))< 0.1],b(t)=1];
slip := [[0, abs(1-x(t))>0.1], b(t) = 0];
sol:=dsolve({sys,b(0)=0},numeric,discrete_variables=[b(t)::boolean],events=[stick,slip]);
plots:-odeplot(sol,[t,x(t)],0..0.01);
plots:-odeplot(sol,[t,x(t)],0..8);
plots:-odeplot(sol,[t,diff(x(t),t)-0.1],0.01..8);


@maple fan This behavior is not a bug.
Take F1.
piecewise first examines if t>=0. If that is the case then it returns 0. It does not go on to the next condition to see if that is also satisfied. It only keeps going if previous conditions are not satisfied.
Thus the plot is a plot of the zero function.
If you want the same result as F2 without using 'and' you could do:
F1a:=piecewise(t<2,0,t<2.1, t^3, exp(t));


@maple fan Here is an answer to the first question.
I use two discrete variables, one for turning on and off and one for getting the time where y(t)=30.
restart;
F1:=y(t)^2;
ode:=diff(y(t),t$2)=F1-b1(t)*piecewise(t<b2(t)+0.1,y(t)^3,exp(y(t)));
ic:=y(0)=1.2,D(y)(0)=0;
res:=dsolve({ode,ic,b1(0)=0,b2(0)=0},numeric,discrete_variables=[b1(t)::boolean,b2(t)::float],events=[[y(t)=30,[b1(t)=1,b2(t)=t]]]);
plots:-odeplot(res,[t,y(t)],0..5);
res(3);
plots:-odeplot(res,[t,b2(t)],0..3);
#############
For the second question I took another example for F1 and introduced a counter b3(t).
There are two events. The first is [y(t)=y0,b3(t)=b3(t)+1], which updates b3(t) when y(t)=y0.
The second is conditional and has a double action: [[y(t)=y0,b3(t)=2],[b1(t)=1,b2(t)=t]].
This latter triggers when y(t)=y0 provided b3(t)=2. If so b1(t) is set to 1 and b2(t) to t.

restart;
F1:=sin(y(t)^2);
ode:=diff(y(t),t$2)=F1-b1(t)*piecewise(t<b2(t)+0.1,y(t)^3,exp(y(t)));
ic:=y(0)=1.2,D(y)(0)=0;
y0:=1.9:
res:=dsolve({ode,ic,b1(0)=0,b2(0)=0,b3(0)=0},numeric,discrete_variables=[b1(t)::boolean,b2(t)::float,b3(t)::integer],
  events=[[y(t)=y0,b3(t)=b3(t)+1],[[y(t)=y0,b3(t)=2],[b1(t)=1,b2(t)=t]]]);
plots:-odeplot(res,[[t,y(t)],[t,y0],[t,y0*b1(t)]],0..8,linestyle=[1,2,1]);
res(4);
plots:-odeplot(res,[t,b3(t)],0..4,thickness=3);




@mahmood180 As I wrote you need the output to be
Array(1..N, 0..M , b)
not what you had (which involved b[m,n]).
Secondly, change the assignment
OB(i,M):=...
to
OB[i,M]:=...

@one_man In copying quotes ``were left out of `minus`.

Converting to set and then to list may change the order:
restart;
a := [y+1, x+2, x+3, x+4]; #Notice the y
a := convert(a, set);
a :=`minus`(a, {x+2});
a := convert(a, list); #Notice the order
###########
## Another example:
restart;
a := [x+5, x+2, x+3, x+4];
a := convert(a, set);
a :=a minus {x+2};
a := convert(a, list);



@sirmotion If you follow the link I gave you will see that 'res' refers to the d'Alembert formula.

@ESchiltkamp I have edited your last worksheet to include the numerical approach.
Notice several points:
1. I use the parameters option for dsolve/numeric, see
?dsolve,interactive,parameters
2. fsolve fails to solve the two equations (really numerical procedures = 0).
But I force it to print out reasonably good solutions.
3. Knowing a reasonably good solution I use that as an initial point for optimizing the sum of squares of the two quantities we want to be zero.
This turns out to give a very good solution.

The worksheet:
MaplePrimes14-10-10odesysAA.mw

Please feel free to ask further questions.

Not answering your actual question just giving the value assuming that the parameters are real and different from zero in which case we may safely assume that they are positive.
A:=Int(exp(-xi^2*b^2)*cos(xi^2*a*t)*cos(xi*x), xi = 0 .. infinity);
value(A) assuming positive;

@ESchiltkamp I'm now very certain that the exact result returned by dsolve for your ode named 'somFy' is wrong. I suspected that, but am now sure. The difference in your case between the correct solution and the dsolve-solution is small enough to make the result look reasonable. However, it is wrong.

It helps to look at a version with simpler constants:
restart;
ode:=-1-diff(y(t), t)^2*signum(diff(y(t), t)) =diff(y(t), t, t);
ics:=y(0) = 0, D(y)(0) = 1;
res:=dsolve({ode,ics});
rr:=odetest(res,ode); #Ought to be zero
plot(rr,t=-Pi/4..3*Pi/4,-.1..1,thickness=3);
#At t = Pi/4 the derivative of res changes sign since cot changes sign at Pi/2. After that rr is nonzero.
diff(res,t);
#The result res is the same as
u:=ln((sin(t)+cos(t))^2)/2;
#thus equal to ln(abs(sin(t)+cos(t))) for t real.
#Now ignoring signum:
ode2:=eval(ode,signum=1);
res2:=dsolve({ode2,ics});
odetest(res2,ode2);
##########
So res and res2 agree when sin(t)+cos(t) > 0. However, that is the case for t=0. Since sin(t)+cos(t) > 0 on -Pi/4 .. 3*Pi/4 and  tends to 0 as t tends to 3*Pi/4 from the left and also as t tends to -Pi/4 from the right the solutions both have maximal domains equal to the open interval -Pi/4 .. 3*Pi/4. That res makes sense for t > 3*Pi/4 is irrelevant because of the singularity at 3*Pi/4.
##########
#Solving ode numerically:
resN:=dsolve({ode,ics},numeric);
evalf(Pi/4);
plots:-odeplot(resN,[t,y(t)],-Pi/4..5); pnum:=%:
##########
#By solving ode2 up to t=Pi/4 (where y ' (t) = 0) and continuing with
ode3:=eval(ode,signum=-1);
#with the approbriate ics at t=Pi/4 it is starightforward to finf the actual exact solution:
RES := y(t) = piecewise(t < Pi/4, ln(sin(t)+cos(t)), t+Pi/4+(3/2)*ln(2)-ln(exp(2*t)+exp((1/2)*Pi)));
##Plotting
plot(rhs(RES),t=-Pi/4..5); pex:=%:
##Displaying RES with the numerical solution
plots:-display(pnum,pex); #They fall on top of each other.

I shall file a bug report (SCR).

@ESchiltkamp Changing alpha to 20 degrees does give us an answer to displacementy which doesn't contain a RootOf. However, I believe it is wrong. First of all it is suspicious that there is no cutoff like what happens in a piecewise. Secondly, odetest gives an identically zero result from t=0 to some t0, after that it continues continuously to nonzero values (untill it blows up, which is not the problem).
displacementy := dsolve({ics2, somFy});
rrr:=odetest(displacementy,somFy): #Testing, rrr ought to be zero
plot(eval(rrr,v0=100),t=0..0.3,0..1);
###############
I tried a different approach. Purely numerical.
Digits:=15:
res:=dsolve(sys,numeric,parameters=[v0],output=listprocedure);
X,Y:=op(subs(res,[sx(t),sy(t)]));
res(parameters=[5]); #Setting the v0 parameter for a test
odeplot(res,[sx(t),sy(t)],0..10,numpoints=100);
odeplot(res,[[t,sx(t)],[t,sy(t)]],0..3,numpoints=100);
odeplot(res,[[t,diff(sx(t),t)],[t,diff(sy(t),t)]],0..3,numpoints=100);
#Now setting up a way to find the speed v0 necessary.
p:=proc(T,v0) option remember; res(parameters=[v0]); [X(T)-distance,Y(T)] end proc;
px:=proc(T,v0) p(_passed)[1] end proc;
py:=proc(T,v0) p(_passed)[2] end proc;
p(1,5); #Testing
sol:=fsolve([px,py],[3,3e12]); #A huge speed guess
res(parameters=sol[2..2]); #For safety I set the parameter v0
odeplot(res,[sx(t),sy(t)],0..sol[1]);
X(sol[1]); #This value ought to be 2.5, but it is somewhat close
#I suspect that there are some numerical problems since the input of [3,3e12] into fsolve returns a listwith the v0-guess unchanged. It is probably due to the huge value of v0.
Trying with
sol:=fsolve([px,py],[3,3.18e12]);
results in a value closer to distance=2.5.
Question: Are the parameters rho and Cd in your problem correct?
You have an enormous resistance requiring an enormous starting speed for this very light ball to hit the ground at 2.5 m.
#################
After seeing your corrected version I tried changing the procedure p above to
p:=proc(T,v0) option remember;
  local L,eps;
  eps:=10^(-Digits);
  res(parameters=[v0]);
  L:=[X(T)-distance,Y(T)];
  if L[1]^2+L[2]^2<eps then print([T,v0],L) end if;
  L
end proc;
px:=proc(T,v0) p(_passed)[1] end proc;
py:=proc(T,v0) p(_passed)[2] end proc;
##That way I get away from the somewhat unwieldy fsolve. I force it to print out results that are good enough for me, but not not good enough for fsolve.
sol:=fsolve([px,py],[0..1,0..7]);


@Carl Love The statement "Anything inside the quotes is treated as a variable, regardless of any other meaning(s) that the characters have." should be modified in view of the results of
`sin`(Pi);
sin(`Pi`);
`sin`:=5; #You won't be allowed to do this
and similar.
The help page ?names says it like this:
"Any valid Maple name formed without using left single quotes is precisely the same as the name formed by surrounding the name with left single quotes. "
sin and Pi are (also) names.

@shiMusa
Yes, the point is not that the numbers should be rationals. The point is that they should not be floats.
The following is straightforward:
restart;
M := Matrix(2, 2, {(1, 1) = a+2.5*I, (1, 2) = 1-I*a, (2, 1) = 4, (2, 2) = a});
M3:=eval(M,a=0.5+0.5*I);
convert(M3,rational);
##
This version is more subtle:
restart;
M := Matrix(2, 2, {(1, 1) = a+2.5*I, (1, 2) = 1-I*a, (2, 1) = 4, (2, 2) = a});
a:=0.5+0.5*I;
M;
convert(M,rational); #No effect
#However, elementwise application of convert works fine:
convert~(M,rational); #Notice the tilde, which implies elementwise application of convert
#To get an idea of what is going on, try debugging:
debug(`convert/rational`);
convert(M,rational);
#Notice the input and the contents of 'flts'.
# Continue with
convert~(M,rational);
#Notice the difference!
#I guess the conclusion is that you should use elementwise conversion.
# To get back to the normal state use undebug:
undebug(`convert/rational`);
##Now try at the element level: first row, first column:
convert(M[1,1],rational);
# We can continue with
indets(M,name);
%;
indets(M,float);



@oldstudent Actually I edited your expression eval( x^3-3*x^2-24*x+8, [x = -2])  to what I wrote in my answer above. Notice in particular the tilde (~) which is used twice in this example. 
The first that happens in
eval~( x^3-3*x^2-24*x+8, x=~[-2,4,6,9]);
is that x=~[-2,4,6,9] becomes [x = -2, x = 4, x = 6, x = 9], which you can try like this:
x=~[-2,4,6,9];
This is use of elementwise `=`.
After that eval is applied elementwise to the list of x-equations as in
eval~(x^3-3*x^2-24*x+8,[x = -2, x = 4, x = 6, x = 9]);


First 138 139 140 141 142 143 144 Last Page 140 of 231