Preben Alsholm

13471 Reputation

22 Badges

20 years, 250 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

[y(t), diff(y(t),t)>0] is a conditional event. The event is that y(t) becomes 0, but under the condition that y'(t) > 0.

Similarly,

events = [[[diff(y(t),t),y(t)<0], halt]]

would make the integration stop when y'(t) = 0, but only if y(t) < 0.

But you are right; the help page for dsolve,events is not easy reading.

You have a typing error. u = t1[i] surely was meant to be u:=t1[i].

I'm wondering (as I was in your previous and related question) why you make substitutions for sin(m*phi) and cos(m*phi). These don't appear in your expression.

Also I think you could simplify the following two commands:

p1:=(beta,Q,P) -> select(x->Im(evalf(eval(x,[beta=beta,Q=Q,P=P])))>0,t1a);
t1:=p1(beta,Q,P):

to just

t1:=select(x->Im(evalf(x))>0,t1a);

Instead of eval(V[x],x=0) =1 you should write

D[1](V)(0,y)=1

which means the partial derivative of V w.r.t. the first variable at (0,y) .

The appearance of a RootOf expression is not necessarily due to a bug. See the help page for RootOf.

You should be able to upload a Maple file by using the icon with a thick green arrow (the last one on the menu).

 

Example.

Try the following two commands.

restart;
solve(d^3-32*(x-l)=0,d);
   
solve({d^3-32*(x-l)=0,t=d},{d,t});
              
In Maple 12 (and also in Maple 15) the first one results in 3 explicit results, whereas the second results in RootOf expressions. This has nothing to do with document mode. In fact I always use worksheet mode.

The presence of the imaginary part after the use of evalf is due to rounding. Try evalf[20](...) .

restart;
s1 := (1/(2*Pi)+Q*cos(k*phi))^2/(1+(1/(2*Pi)+Q*cos(k*phi))^2);

Here is my version built upon yours.
Q:=15;
#The integral from 0 to 2*Pi is what you want and that is independent of the integer k, so:

k:=1;
#The integral over 0..2*Pi is twice the integral over 0..Pi
A:=Int(s1,phi=0..Pi);
evalf(A);
#Changing variables.
IntegrationTools:-Change(A,phi=2*arctan(u),[u]);
evalindets(%,specfunc(anything,cos),simplify@expand);
normal(%);
s2:=IntegrationTools:-GetIntegrand(%);
t:=singular(s2,u);
t1a:=map(rhs@op,[t]);
t1:=select(x->evalb(evalf(evalc(Im(x))>0)),t1a);
Res:=add(residue(s2,u=p),p=t1):
simplify(2*Pi*I*Res);
evalf[20](%);

If on the surface of it you see common factors you could force factorization by a procedure like the following.

sfactor:=proc(u) local i,F,S,p;
if not hastype(u,`+`) then return u end if;
if not type(u,`+`) then return evalindets(u,`+`,procname) end if;
for i to nops(u) do
  if not type(op(i,u),`*`) then
    F[i]:={op(i,u)}
  else
    F[i]:={op(op(i,u))}
  end if
end do;
S:=`intersect`(seq(F[i],i=1..nops(u)));
if S={} then
   u
else
   p:=`*`(op(S));
   p*map(`/`,u,p)
end if
end proc;

#Test on your example.
u:=(Pi^2*(c+(1+Pi)/Pi)^2+Pi^2)/((c-(-1+Pi)/Pi)^2*Pi^2+Pi^2);
T:=op(1,u);
sfactor(T);
sfactor(u);

What takes time is evaluating the sum. When k = 0 that time is wasted, since the result in your particular example obviously is zero.

In your example you could have proceeded in one of the following two ways.

#Method 1: unevaluation quotes around s.
restart;
f:=(k,y)->k*'s'(y);
s:=y->add(sin(1/y^i),i=1..1000000);
f(0,y);
f(1,y);
%;
#Method 2: using the inert form for add, %add.
restart;
f:=(k,y)->k*s(y);
s:=y->%add(sin(1/y^i),i=1..1000000);
f(0,y);
f(1,y);
value(%);

With no details about the more general cases your mention it is impossible to say if these ideas are useful.

In the help page for SingularValues you find:

"The singular values S together with the left singular Vectors as columns of Matrix U and the right singular Vectors as columns of Matrix Vt satisfy
                       (U . SS) . Vt = A
. Matrix SS has the entries of S along its main diagonal."

Example:
restart;
with(LinearAlgebra):
A := RandomMatrix(5,3,datatype=float);
S,U,Vt:=SingularValues(A,output=['S','U','Vt']);
S1:=DiagonalMatrix(S[1..3],5,3);
U.S1.Vt;
Equal(evalf(%),A);

If you are not already aware of dualaxisplot from the plots package, you may want to explore the possibilities it offers together with options for label directions and gridlines.

You can try removing the initial condition:

remove(has,edeq,_C[0]);
dsolve(%);

You will see (this is Maple 12 at the moment) that Maple cannot find an explicit expression for the general solution, so it is no wonder that solve must return NULL when attempting to solve the equation involved in satisfying the initial condition.

I assume that you meant

A:= subs(Fsubs,parvals,Theta0/C=(exp(Pec_i*F*(1-S0))-1)*(H+1/(1-exp(-C/k*Pec_i*F*S0))));

B is not an equation and seems to be missing a parenthesis. If Maple expects an equation and don't get one, often it assumes you meant equal to zero. Probably you didn't mean that.

With the modification of A above you have a connection between R3  and S0. You can explore that by doing

plots:-implicitplot(A,R3=-3..3,S0=-2..2,gridrefine=2);

Suppose you are interested in the case S0 >0, then you could define numerically S0 as a function of R3 like this:

f:=r->fsolve(eval(A,R3=r),S0=1);

#Test and plot
f(0.12345);
plot(f,0..3);


As long as you have an initial value problem you shouldn't have much of a problem. With the conditions at tf you have an initial value problem.

I have access to Maple 12 right now so will make use of map instead of the elementwise operation ~ which is available in Maple 13 and later.

Using MatrixExponential from the LinearAlgebra package:

with(LinearAlgebra);
A:=Matrix([[-1/R/C,1/C,0,0],[-1/L,0,0,0],[-2,0,1/R/C,1/L],[0,0,-1/C,0]]);
b:=<0,0,2*nf,0>;
Eigenvalues(A);
W:=unapply(MatrixExponential(A,t),t):
X:=W(t-tf).<x1f,x2f,p1f,p2f>+map(int,W(t-s).b,s=tf..t);

Using dsolve (and making use of A and b above.

xp:=<x1(t),x2(t),p1(t),p2(t)>;
S:=map(diff,xp,t)=A.xp+b;
sys:=convert(lhs(S)-rhs(S),set);
#dsolve without specifying the method gives a messy result:
res1:=dsolve(sys union {x1(tf)=x1f,x2(tf)=x2f,p1(tf)=p1f,p2(tf)=p2f}):
# method=laplace is quicker and neater.
res2:=dsolve(sys union {x1(tf)=x1f,x2(tf)=x2f,p1(tf)=p1f,p2(tf)=p2f},{x1(t),x2(t),p1(t),p2(t)},method=laplace);

In Maple 15 the following works.

restart;
sys:={diff(b1(t), t) = I*u1(t),
diff(b2(t), t) = I*u2(t)-3*b1(t)/(2*omega),
diff(b3(t), t) = I*u3(t),
diff(u1(t), t) = I*b1(t)+2*u2(t)/omega-gradQ1,
diff(u2(t),t)=I*b2(t)-1/(2*omega)*u1(t) - gradQ2,
diff(u3(t), t) = I*b3(t)+I*b3(t)-gradQ3};
gradQ2:= (3*K*u2(t)-2*u1(t)*omega)/((9/4)*K+omega^2+omega^2/kappa^2);
gradQ1 := 1.5 * K * gradQ2/omega;
gradQ3 := gradQ2/kappa;
K := (x + t);
inits:={b1(0) = exp(-(1/25)*ln(2)*x^2), b2(0) = (3/25*I/omega)*ln(2)*x*exp(-(1/25)*ln(2)*x^2), b3(0) = 0, u1(0) = 0, u2(0) = 0, u3(0) = 0};
res:=dsolve(sys union inits,numeric,parameters=[x,kappa,omega],range=0..16,output=listprocedure);
B1,B2,B3,U1,U2,U3:=op(subs(res,[b1(t),b2(t),b3(t),u1(t),u2(t),u3(t)]));
#Test
res(parameters=[1,.125,.125]);
B1(2);
#Now we make a procedure accepting values for the 3 parameters and for t.
#It has been made to return unevaluated if the parameters are not of type numeric.
#B1(t) itself will return unevaluated if t is not of type numeric.
Qb1:=proc(xx,k,w,t) if not type([xx,k,w],list(numeric)) then return 'procname(_passed)' end if;
   res(parameters=[xx,k,w]);
   B1(t)
end proc;
#Test
Qb1(x,.125,.125,t);
Qb1(1,.125,.125,t);
Qb1(1,.125,.125,2);
with(plots):
complexplot(Qb1(1,.125,.125,t),t=0..1);
animate(complexplot,[Qb1(x,.125,.125,t),t=0..0.5,thickness=3],x=-1..1,frames=10);

#When you assigned 'res' above you may have noticed that Maple didn't print anything to the screen.
This seems like a bug to me since it will normally do that. What exactly triggers this I don't know, but have noticed that the the following command, where I is  simply replaced by 1, prints the expected.

Res:=dsolve(eval(sys union inits,I=1),numeric,parameters=[x,kappa,omega],range=0..16,output=listprocedure);


#I also tried the whole code in Maple 12 and it didn't work, but that is not surprising since even the parametrized initial value problem (notice the presence of exp(a), not just a or e.g. a^2).

eq:=diff(y(t),t)=y(t):
res2:=dsolve({eq,y(0)=exp(a)},numeric,parameters=[a]);

#gives the error
Error, (in dsolve/numeric) invalid specification of initial conditions

Added. I have submitted a bug report (aka "Software change request"). The printing problem has nothing to do with using the parameters option, but seems to occur because of the presence of imaginary numbers in the ivp.

Here is a solution to the first system using the procedure option in dsolve/numeric
?dsolve,numeric,ivp

restart;
odes := diff(x(t), t, t)*sin(diff(y(t), t))+diff(y(t), t, t)^2 = -2*x(t)*y(t)*exp(-diff(x(t), t))+x(t)*diff(y(t), t)*diff(x(t), t, t), x(t)*diff(y(t), t, t)*diff(x(t), t, t)+cos(diff(y(t), t, t)) = 3*y(t)*diff(x(t), t)*exp(-x(t)):
ics := x(0) = 1, (D(x))(0) = 0, y(0) = 0, (D(y))(0) = 1:
#We rewrite the system as a first order system:
subs(diff(x(t),t,t)=yp2,diff(y(t),t,t)=yp4,diff(x(t),t)=Y[2],diff(y(t),t)=Y[4],x(t)=Y[1],y(t)=Y[3],{odes});
#The result of this is pasted as the equations 'eqs' in the procedure.
#YP[2] and YP[4] are found by fsolve.
p:=proc(N,t,Y,YP) local eqs,yp2,yp4;
   YP[1]:=Y[2];
   YP[3]:=Y[4];
   eqs:={yp2*sin(Y[4])+yp4^2 = -2*Y[1]*Y[3]*exp(-Y[2])+Y[1]*Y[4]*yp2,
        Y[1]*yp4*yp2+cos(yp4) = 3*Y[3]*Y[2]*exp(-Y[1])};
   YP[2],YP[4]:=op(subs(fsolve(eqs,{yp2,yp4}),[yp2,yp4]));
end proc:
res:=dsolve(numeric,procedure=p,initial=Array([1,0,0,1]),number=4,procvars=[x(t),diff(x(t),t),y(t),diff(y(t),t)],start=0);
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..5);

Inside procedures you don't have full evaluation unlike the situation for commands from the input line.

That could be the problem, but if we see your definition of `*` we may be able to settle that question.

First 131 132 133 134 135 136 137 Last Page 133 of 158