Preben Alsholm

12103 Reputation

22 Badges

16 years, 196 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

I tried your worksheet in Maple 18, Maple 2020.2, and in Maple 2021.1.

The double loop over i and j got stuck in Maple 2020.2, but not in the other two.

Already the computation of evalf(int(r*h0[i]*g[j], r = 0 .. 1)) gets stuck for i = j = 1.
Try this version instead:
 

for i to n do
    for j to n do h[j][i] := int(r*h0[i]*g[j], r = 0 .. 1,numeric,method=_d01ajc)/int(r*(-r^2 + 1)*g[j]^2, r = 0 .. 1,numeric,method=_d01ajc); end do;
end do;

The use of method=_d01ajc makes the computation real fast.

Try this:
 

expr:=piecewise(x(t)<0, -1, x(t)<1, 0, 1) + f(t) + piecewise(t<0, 0, t<1, 1, 0) + x(t)*piecewise(t<0, 1, t<1, 0, 1) + a*piecewise(t<0, 3, t<1, -1, 2);
####
select(hastype,expr,specfunc(piecewise));
select(hastype,%,And(relation,satisfies(r->lhs(r)=t)));

 

Trying this:

plot(w->int(eval(A1+A2,omega=w),theta=0..2*Pi,numeric,digits=5),0..50);

I got

You won't be able to find a symbolic solution. You must solve numerically.

Thus you must give values to the parameters c and t.
You must also give boundary or intial values.

I'm assuming that by e^() you mean the exponential function, which in Maple is exp. So I have made that replacement.
Here is an example of solving the initial value problem.

restart;
sys := -diff(g(x), x)*exp(c*x)*c/(3*t) + f(x) = 0, t*(c*exp(2*c*x)*diff(g(x), x) + t*exp(c*x)*(2*c^2 + 3*f(x) - 1))*diff(f(x), x) + 6*c^2*(diff(f(x), x, x)*t*exp(2*c*x) + diff(f(x), x, x, x)*exp(3*c*x)/3) = 0;
sol:=dsolve({sys}); # essentially unsolved, but work has been done
dsolve(sol[1]); # returns NULL i.e. nothing
#Solving the original system numerically:
res:=dsolve({sys,f(0)=f0,D(f)(0)=f1,(D@@2)(f)(0)=f2,g(0)=g0},numeric,parameters=[c,t,f0,f1,f2,g0]);
#Setting parameters (an example):
res(parameters=[1,2,0.5,0.6,0,0]);
plots:-odeplot(res,[[x,f(x)],[x,g(x)]],-1..10);
## Here is a procedure that can be used by animate:
Q:=proc(c,t,f0,f1,f2,g0,scene::list:=[[x,f(x)],[x,g(x)]],x_interval::range:=-1..10) 
   if not [c,t,f0,f1,f2,g0]::list(realcons) then return 'procname(_passed)' end if;
   res(parameters=[c,t,f0,f1,f2,g0]);
   plots:-odeplot(res,scene,x_interval,_rest)
end proc;
## Examples of use all with c, t, f1, f2, g0 fixed
plots:-animate(Q,[1,2,f0,0.6,0.1,-.3,0..20],f0=-1..1,trace=24); 
plots:-animate(Q,[1,2,f0,0.6,0.1,-0.3,[f(x),diff(f(x),x)],0..30],f0=-1..1);
plots:-animate(Q,[1,2,f0,0.6,0.1,-0.3,[seq([x,diff(f(x),[x$k])],k=0..2)],-0.3..5],f0=0..5);

In the worksheet attached I end up displaying all the separatrices in sequence producing the animation below.
MaplePrimes21-05-05_separatrices.mw

Since the problem is formulated in vector form the Maple solution should be given as such, I think.
Even Maple 2021 doesn't quite do that.
Here is the closest I got:
 

restart;
X:=<x1(t),x2(t)>;
Sys:=diff(X,t)=<-4,-2;6,3>.X+<2/(exp(t)+1),-3/(exp(t)+1)>; 
ics:=eval(X,t=0) = <0,ln(2)>;
 
sol1:= dsolve([ Sys, ics]); # Not on vector form
sol:=X=eval(X,sol1); # Now on vector form

The solution looks like this:

You have a classical Volterra-Lotka system.
All solutions are periodic. The orbits in phase space are closed.
Use odeplot with a numerical solution as done here.
 

restart;
sys := diff(x1(t), t) = 3*x1(t) - 2/1000*x1(t)*x2(t), diff(x2(t), t) = 6/10000*x1(t)*x2(t) - 1/2*x2(t);
ics := x1(0) = 1000, x2(0) = 500;
res:=dsolve({sys,ics},numeric);
plots:-odeplot(res,[[t,x1(t)],[t,x2(t)]],0..15,legend=[x1,x2],size=[1000,500]);
plots:-odeplot(res,[x1(t),x2(t)],0..7,frames=100);

The last is an animation of the orbit in phse space and is seen below.

I don't have Maple 13, but I have Maple 12. It gives the same error.
The error is NOT related to the way you produce the subscript, which is seen by this code:
 

restart;
H1:=v*P[r]-((1-P[r]^2*r/(P[r]^2*r-1))/r)^(1/2);
simplify(H1); ## Error in Maple 12
normal(H1);   ## OK
radnormal(H1);         ## Same error
simplify(H1,symbolic); ## Same error

The error is also present in Maple 15, 16, 17, 18, and 2015, but not in Maple 2016.

@nm For one session you can can keep the order once determined by giving LinearAlgebra:-Eigenvectors option remember or maybe option remember,system:
 

LA_EV:=subsop(3=(remember,system),eval(LinearAlgebra:-Eigenvectors)):

So you could compare the outputs of LA_EV and LinearAlgebraEigenvectors like this:
 

restart;
Sy:=1/sqrt(2)*Matrix([[0,-I,0],[I,0,-I],[0,I,0]]);
LA_EV:=subsop(3=(remember,system),eval(LinearAlgebra:-Eigenvectors)):
to 5 do LA_EV(Sy) end do;
to 5 do LinearAlgebra:-Eigenvectors(Sy) end do;

This only assures that if you happen to execute another time the previous result will be returned.
Between sessions this is no solution. You must sort.
 

Replace expand by solve or fsolve.
Here I use fsolve. Notice that all 3 roots are real (as it happens):
 

restart;
A1 := x^3 + n - 81*x - 6 = 0;
for n to 11 do
    fsolve(A1)
end do;

 

Since Z is taking values in the complex plane and E is therefore also, you would need a 4D coordinate system to illustrate E's dendence on Z. Maple uses color as the fourth dimension in complexplot3d. Take a look at the help page.

rho := 1/4; mu := 1/4; E := add(Z^k/GAMMA(k*rho+mu), k = 0 .. 5);
plots:-complexplot3d(E,Z=-5-5*I..5+5*I);

You simply need to give epsilon a numeric value. Then your code runs.
epsilon is present in exz, which is going to be printed.

@AHSAN In your worksheet solve gives its answer as a RootOf expression.
That expression is not at all useless.
 

restart;
pg1 := -(888888889*sigma^6*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(5/2) - 888888889*sigma^6*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(5/2) + 3333333333*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*k*sigma^5 + 3333333333*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*k*sigma^5 + 6666666666*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*Q*sigma^4 - 3333333333*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*sigma^5 + 6666666666*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*Q*sigma^4 - 3333333333*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*sigma^5 + 15000000000*k^2*lambda*sigma^2 + 60000000000*Q*k*lambda*sigma - 30000000000*k*lambda*sigma^2 + 60000000000*Q^2*lambda - 60000000000*Q*lambda*sigma + 15000000000*lambda*sigma^2)/(5000000000*(k*sigma + 2*Q - sigma)^2*sigma^3);
S := solve(pg1 = 0, sigma);
Sa:=allvalues(S):
numelems({Sa});
evalf(eval([Sa],{k = 2, Q = 3, lambda = 4}));
select(type,%%,realcons);

 

While evalc assumes that any variables present are real, convert/trig doesn't.
In fact the equation exp(I*z) = cos(z) + I*sin(z) holds for all z in the complex plane.

evalc(exp(I*x));
convert(exp(I*x),trig);

 

Here you will see the solution method:
 

restart;

ode:=diff(y(x),x)*tan(diff(y(x),x))+ln(cos(diff(y(x),x))) = y(x);
infolevel[dsolve]:=5:
dsolve(ode);
DEtools:-odeadvisor(ode); # Wrong

The printout is:

Methods for first order ODEs:
-> Solving 1st order ODE of high degree, 1st attempt
trying 1st order WeierstrassP solution for high degree ODE
trying 1st order WeierstrassPPrime solution for high degree ODE
trying 1st order JacobiSN solution for high degree ODE
trying 1st order ODE linearizable_by_differentiation
trying differential order: 1; missing variables
<- differential order: 1; missing "x" successful

 

1 2 3 4 5 6 7 Last Page 1 of 151