Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@BrettKnoss You can do that in this way:
 

restart;

f:= 2/(3-x);

s:=convert(f,FormalPowerSeries);

term__k:=simplify(op(1,s)); 

subs(x=abs(x),term__k);
ratio:=subs(k=k+1,%)/%;
limit(ratio,k=infinity);  # answer abs(x)/3

Thus the radius of convergence is 3.

@nm I tried your original worksheet again in Maple 2020.2, Maple 2019.2, and Maple 2015.2.
The first two gave the simple answer twice in two runs.

Maple 2015.2 gave in two runs: (1) simple, not so simple and (2) not so simple, not so simple.

@nm I updated my Physics package.
Still no difference between warnlevel = 3 and warnlevel = 0.

@gasmiboubeker I cannot reproduce the problem in Maple 18.02 nor in Maple 2018.2.

I used the code as written by Tom Leslie.

My Physics version at the time of my first (and so far only) run on your worksheet is 879, whereas yours is the most recent i.e. 884.
With my version the result looked exactly the same in the two cases. The answer looked like the second of yours, i.e. with warnlevel = 3 (the default).
I ran the worksheet again with the same result.
I hesitate to update my Physics because of this. At least I shall wait a few hours.

`Standard Worksheet Interface, Maple 2020.2, Windows 10, November 11 2020 Build ID 1502365`

`The "Physics Updates" version in the MapleCloud is 884. The version installed in this computer is 879 created 2020, November 19, 20:10 hours Pacific Time, found in the directory C:\\Users\\Bruger\\maple\\toolbox\\2020\\Physics Updates\\lib\\`

@Carl Love When I copied your code I got a syntax error. There are more left parentheses than right.
So it should be:
 

ode2:= -v+mu*diff(r*diff(u(r), r), r)/r = 0;

 

@nm The statement in the help that you quote:

"In some cases, the execution may not abort at exactly the time limit imposed, but will abort as soon as it can do so safely"

seems to be contradicted in this example by the fact that the user can stop the computation by pressing the red icon on top.
That is not always the case though.

Why is your title :"Maple write differential equation from Lotka-Volterra Model" ?

This has absolutely nothing to do with Lotka-Volterra.

@mmcdara Classical methods in dsolve (like method=classical[rk4]) use fixed stepsize and has no error control. That allows us to use your original piecewise version mentioned in your original question.
Results appear good with stepsize = 1e-4.
 

# Using P_DATA100 the system is: 
sys100 := {10*diff(v(t), t) = 900*t - 1000 - 10000*x(t) - 100*piecewise(0 < v(t), 1, v(t) < 0, -1, 0), diff(x(t), t) = v(t), v(10/9) = 0, x(10/9) = 0};
#Now using Runge-Kutta 4:
res:=dsolve(sys100,numeric,method=classical[rk4],stepsize=0.0001);
odeplot(res,[t,x(t)],T__0..2);

@mmcdara I have attached a successful events version.
I will say, however, that the classical approach is simpler and real fast with stiff=true (i.e. rosenbrock).
The discrete variable f(t) used here is of type float and will take on 3 values: -1, 0, or 1.
Events are triggered by v(t) +eps =0, v(t)-eps = 0, and v(t)=0.
For both methods stiff=true is chosen although it probably has no effect in the events case.

 

restart;

with(plots):

# x(t)        body displacement
# v(t)        body velocity
# M__p     body mass
# V__p     body volume
# P__0      pre-stress of the spring
# K           spring stiffness
# C           solid-solid friction coefficient
# rho__h   density of the fluid the body is immersed in
#
# gamma__longi(t) longitudinal acceleration

gamma__longi := t -> A*t;

####################################

# The classical approach:

equationsC := {
               diff(x(t), t) = v(t),
               M__p * diff(v(t), t) = (M__p - V__p*rho__h) * gamma__longi(t)
                                      - (P__0 + K*x(t))
                                      - C*tanh(v(t)/v__0)
            };

 

EVENTS:=[[v(t)+eps,f(t)=-1],[v(t)-eps,f(t)=1],[v(t),f(t)=0]];

equationsE := {
               diff(x(t), t) = v(t),
               M__p * diff(v(t), t) = (M__p - V__p*rho__h) * gamma__longi(t)
                                      - (P__0 + K*x(t))
                                      -f(t)*C
            };

# take-off time
# (obtained by equating to 0 the rhs of the 2nd edo after having set x(t)=v(t)=0 and removed the friction term)


rest_state := eval(equationsC[1], [x(t)=0, v(t)=0]);
T__off     := solve(eval(rest_state, C=0), t);
 

# initial conditions

IC_C := { x(T__off) = 0, v(T__off) = 0};
IC_E := { x(T__off) = 0, v(T__off) = 0, f(T__off) = 0};
# full systems

sysC := equationsC union IC_C;
sysE := equationsE union IC_E;

eps:=1e-4: v__0:=1e-7:

SOL_C:=dsolve(sysC,numeric,parameters=[A, M__p, V__p, rho__h, P__0, K, C],stiff=true); ## stiff
SOL_E:=dsolve(sysE,numeric,parameters=[A, M__p, V__p, rho__h, P__0, K, C],
                   discrete_variables=[f(t)::float],
                   events=EVENTS,maxfun=10^6, stiff=true  #stiff doesn't help much here
                );

P_DATA0 := [
            A      = 100,
            M__p   = 10,
            V__p   = 1,
            rho__h = 1,
            P__0   = 1000,
            K      = 10000,
            C      = 0
          ]:

T__0 := eval(T__off, P_DATA0);
 

SOL_C(parameters=P_DATA0):
Frictionless_C := odeplot(SOL_C, [t, x(t)], t=T__0..2, color=blue);
 

SOL_E(parameters=P_DATA0):
Frictionless_E := odeplot(SOL_E, [t, x(t)], t=T__0..2, color=blue);

P_DATA100 := [
            A      = 100,
            M__p   = 10,
            V__p   = 1,
            rho__h = 1,
            P__0   = 1000,
            K      = 10000,
            C      = 100
          ]:

T__0 := eval(T__off, P_DATA100);
 

SOL_C(parameters=P_DATA100):
t0:=time():
StrongFriction_C := odeplot(SOL_C, [t, x(t)], t=T__0..2, color=red):
time()-t0; # 0.015

display(Frictionless_C, StrongFriction_C); pC:=%:

SOL_E(parameters=P_DATA100):
t0:=time():
StrongFriction_E := odeplot(SOL_E, [t, x(t)], t=T__0..2, color=red):
time()-t0; # 2.250s

display(Frictionless_E, StrongFriction_E); pE:=%:

display(pC,pE);

 


pC and pE look the same and overlap totally in the last display with eps=1e-4.
With eps = 1e-3 computation time is much smaller and overlap is still good.
With eps = 1e-2 computation time is even smaller, but there is clear visual disagreement between pC and pE.
odeplot(SOL_E,[t,f(t)], T__0..2,numpoints=10000) will show you that an enormous amount of events takes place, so no wonder it takes time.

Download MaplePrimes20-11-07_classical_and_event.mw

@Muhammad Usman Commenting out what may interfere with the timing (such as remembering values from the plots) and omitting printing the loop I get 0.187s.
fnum can be seen afterwards simply by eval(fnum);
 

restart;

U := q^6*sin(p);
alpha := 1/3;
M1 := 3;
M2 := 3;
A:=Int(Int(U/((x^rho - p^rho)^alpha*(y^rho - q^rho)^alpha), p = 0 .. x), q = 0 .. y);

Q:=proc(r,xx,yy) evalf(eval(A,{'rho'=r,x=xx,y=yy})) end proc;
##plot3d(Q(1,x,y),x=0..1,y=0..1); # Test
## An animation in rho with only 6 frames
##plots:-animate(plot3d,[Q(rho,x,y),x=0..1,y=0..1],rho=1..6,frames=6);

## The revised loop:
t0:=time():
##printlevel:=3:
for i to 6 do
    for i11 to M1 do
        for j11 to M2 do 
            fnum[i, i11, j11] := Q(i,i11/M1,j11/M2)
        end do
    end do
end do:
time()-t0;
# time = 0.187s
eval(fnum);

 

@mmcdara I just tried these events:
 

events=[[[0,v(t)>0],b(t)=1],[[0,v(t)<0],b(t)=0]]

They seem to produce exactly what you have in your event version.

@mmcdara I would like to, but haven't seen a way out yet.

When I do ?ConvertIn the help page for modp1 opens.
On that page it says:
"The modp1 function accepts the following functions, whose names are also protected global names."

Then follows 5 colums of names some of which are in black (i.e. no link to help page). Among those is ConvertIn.
Its meaning is understood by modp1 and is described on that help page. In fact ConvertIn may just be nothing but a protected name used in a call to modp1 to tell it what to do.
Some of those with an actual link (the blue ones) like Add just brings you to the help page for add. That must be a mistake.
How Add is used is beyond me.

 

A common way of solving this kind of ode is to multiply by diff(y(x),x) on both sides and integrate like this:
 

restart;
ode:=2*diff(y(x),x$2)=exp(y(x));
IC:=y(0)=0,D(y)(0)=1;
diff(y(x),x)*ode;
ode1:=map(int,%,x)+(0=C1);
eval[recurse](convert(ode1,D),{x=0,IC});
## So C1 = 0.
ode2:=eval(ode1,C1=0);
odes:=solve(ode2,{diff(y(x),x)});
ode3:=op(odes[1]);
sol1:=dsolve(ode3);
eval[recurse](sol1,{x=0,IC});
## so _C1 = -2.
sol:=eval(sol1,_C1=-2);

 

First 22 23 24 25 26 27 28 Last Page 24 of 229