Preben Alsholm

12813 Reputation

22 Badges

17 years, 304 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@nm You are right. Something strange is going on.
If in a new worksheet I start with my integrand2 and then do int(integrand2,x) I get an error.
If I then execute simplify(integrand2) and then int(integrand2,x) I get the result.
And if again in a fresh worksheet I start with your original integrand, then do simplify(integrand) then finally
int(integrand,x) then I get a result. Notice that I'm not using the simplified version:

restart;
integrand:=(((-3*x^2-18*x-27)*exp(2)^2+(30*x^3+330*x^2+1170*x+1350)*exp(2)-75*x^4-1200*x^3-7050*x^2-18000*x-16875)*ln(x)+(12*x^2+54*x+81)*exp(2)^2+(-120*x^3-1106*x^2-3510*x-4050)*exp(2)+225*x^4+3560*x^3+20990*x^2+54000*x+50625)/((3*x^4+18*x^3+27*x^2)*exp(2)^2+(-30*x^5-330*x^4-1170*x^3-1350*x^2)*exp(2)+75*x^6+1200*x^5+7050*x^4+18000*x^3+16875*x^2);
simplify(integrand);
int(integrand,x);

I get the errors also from your Calendar version though.
Using restart doesn't remove everything (see ?restart).
So in testing this we must start with a fresh worksheet (a new engine).

 

@ You are wrong. The general term tends to 0:
 

limit(1/k^(2-cos(1/k)),k=infinity); # 0
## Obvious also since the term is positive and satisfies:
1/k^(2-cos(1/k))<= 1/k; 

 

@Kitonum If the result isn't obvious it may be a good idea to see how Maple finds that:
 

restart;
infolevel[int]:=5;
int(1/k^(2-cos(1/k)), k=1..infinity);
# So try this
asympt(1/k^(2-cos(1/k)),k);

So the result infinity appears to be OK.

@dearcia You must have an older Maple version.
Replace DEtools:-DEplot by DEtools[DEplot] and you should be OK.
If, however, your version is older than Maple 13 then elementwise application (using the tilde ~) is not available. You can work around that rather easily by using map.
###############################
Here is a version that I made for Maple 12 on the basis of my Maple 2021 version.
That made me realize how convenient the elementwise operation is. It came with Maple 13.
There are several changes as you will notice.
## This version also works in Maple 2021.
MaplePrimes21-10-22_ode_phaseplot_VersionM12.mw
 

restart;

F0 := 0; zeta := .25; w := 1; Omega := 1; m := 1;

ode1 := diff(x(t), t) = y(t); ode2 := diff(y(t), t) = F0*cos(Omega*t)/m-2*zeta*w*y(t)-w^2*x(t);

res:=dsolve({ode1, ode2,x(0)=x0,y(0)=y0},numeric,parameters=[x0,y0]);

P:=proc(xy0::list,tR::range:=0..10) if not xy0::list(realcons) then return 'procname(_passed)' end if;
   res(parameters=xy0);
   plots:-odeplot(res,[x(t),y(t)],tR,_rest)
end proc;

#Test of P

P([0,50],linestyle=dash,color=blue,thickness=3);

P([0,50],-5..15,linestyle=dash,color=blue,thickness=3);

L:=[[x(0) = 0, y(0) = 50], [x(0) = 9, y(0) = 25], [x(0) = 85, y(0) = 20], [x(0) = .25, y(0) = .5], [x(0) = 7, y(0) = 5]];

L1:=map(x->map(rhs,x),L);

opts:=zip(`=`,[color$5],[red, green, black, navy, maroon]),[(thickness=3)$5],zip(`=`,[linestyle$5],[dash, dashdot, dot, longdash, spacedash]);

nopts:=nops([opts]);

p0:=DEtools[DEplot]([ode1, ode2], [x(t), y(t)], t = 0 .. 1, x = -100 .. 100, y = -100 .. 100,
 [[x(0) = 0, y(0) = 0]],arrows = medium,color=green):

plots[display](seq(P(L1[i],seq(opts[j,i],j=1..nopts)),i=1..5),p0);

plots[display](seq(P(L1[i],-0.9..10,seq(opts[j,i],j=1..nopts)),i=1..5),p0);

 


 

Download MaplePrimes21-10-22_ode_phaseplot_VersionM12.mw

 

@tomleslie OK, I misread the result from A3 as having a minus Float(infinity).

But my real mistake was to use a wrong ordering of the arguments of f3.
I should have considered
 

W1:=eval(W,[m,s,t,w]=~[1, 3, -2, -1]);
plot(Re(W1),x=0..1);

instead of

W1:=eval(W,[s,t,w,m]=~[1, 3, -2, -1]);

which gives no problems.
The plot of the correct W1 shows a clear indication of a singularity at x = 1/3
 

singular(W1,0..1); # {x = 0}, {x = 1/3}, {x = 1}
evalf(Int(Re(W1),x=0..1/3)); # -0.002885966034
evalf(Int(Re(W1),x=1/3..1)); # Problems

 

@tomleslie Obviously the real value is wrong since it must be larger than -2 if the following plots can be believed:
 

int(-1/(1-x)^(1/2),x=0..1);
plot([Re(eval(W,{s=1,t=3,w=-2,m=-1})),-1/(1-x)^(1/2)],x=0..1);

 

@tomleslie laurent only exists as numapprox:-laurent and as such you can see that it relies entirely on the builtin procedure series.

showstat(numapprox:-laurent);

'laurent' is also a type:

showstat(`type/laurent`);

@pazduha Your k can be written simpler:
 

k2:=piecewise(t < t1, kmax, t < 2*t1 + sqrt(2)*t1, -kmax, t < 3*t1 + 2*sqrt(2)*t1, kmax, t < 4*t1 + 2*sqrt(2)*t1, -kmax, 0);

piecewise works this way: If the first condition is satisfied then you get the value given in the second argument no matter what the other conditions are.
If the first condition is not satisfied, then the next condition is considered: If true then the value following that condition is returned, if not then on to the third condition etc.
Since your times t1, 2*t1 + sqrt(2)*t1, 3*t1 + 2*sqrt(2)*t1, 4*t1 + 2*sqrt(2)*t1 are increasing (when t1 >0) then the "And" you are using is unnecessary.
Finally if the last condition (here  t < 4*t1 + 2*sqrt(2)*t1) is not satisfied, you don't need an extra condition together with a value. Just give the value (here 0). You can if that value is actually 0 leave it out; it will be assumed to be zero.

@Muhammad Usman
(1) Your triple loop begins with m = n = k = 0. thus exp(k*eta[3] + m*eta[1] + n*eta[2]) = 1.
     this explains the error.
(2) The first argument to coeff is supposed to be a polynomial in a variable or expression (your case).
     Your Eq1 is a fraction: try numer(Eq1) and denom(Eq1).
     So Eq1 is not a polynomial in exp(k*eta[3] + m*eta[1] + n*eta[2]) for any m,n,k with m+n+k>0.

You should let us know what you are actually looking for.

Eq2 is not defined in that worksheet.

@rcorless Actually odetest also works:
 

odetest(y(t)=theta,diff(y(t),t,t)+omega^2*sin(y(t))=0);  # 0

where theta is defined as in your worksheet.

@rcorless I managed to get zero by treating the two terms in the ode differently:
 

restart;
k := sin(alpha/2);
theta := 2*arcsin(k*JacobiSN(omega*t, k));
t1:=expand(omega^2*sin(theta));
t2:=simplify(diff(theta, t, t));
simplify(t1+t2);  # 0

 

@Muhammad Usman XJTU Since you say that you are new to Maple and have little knowledge of it, you probably shouldn't have just taken over the code written by somebody else for another ode.
Your ode as presented in Maple code is first order, so as pointed out by Tom Leslie and Rouben Rostamian you are only allowed one condition and shooting doesn't make any sense.
Rouben suggests that you present your problem in words and not in Maple code.
I think you should do that. Is your intended ode actually the one we see?

You don't seem to be defining the function PHIb anywhere. What is it?

Perhaps needless to say, but the user needs kernelopts(floatPi=true) as the code is written. That setting is the default.
I have chosen not to use it though,so I just replace 2.*Pi by twoPi, where twoPi is
local twoPi:=evalf(2*Pi).

3 4 5 6 7 8 9 Last Page 5 of 214