vv

6117 Reputation

18 Badges

3 years, 327 days

MaplePrimes Activity


These are replies submitted by vv

@Rouben Rostamian  
Yes, it works, not very intuitive (without a plot) though.

restart;
a:=1/3: h:=1/20: d:=1/4:
M:= # Moebius
<(1+u*a * cos(t/2))*cos(t), (1+u*a * cos(t/2))*sin(t), a*u * sin(t/2)>:
M_u, M_t := <diff(M,u)>, <diff(M,t)>:
NN:=LinearAlgebra:-`&x`(M_u,M_t):
N:=LinearAlgebra:-Normalize(NN,2):
p1:=plot3d( M+h*N, t=0..2*Pi, u=-1..1, scaling=constrained, color=orange):
p2:=plot3d( M-h*N, t=0..2*Pi, u=-1..1, scaling=constrained, color=green):
narr:=10:
q1:=seq( plots:-arrow(eval(M+h*N, [u=0,t=k*Pi/narr]),eval(d*N, [u=0,t=k*Pi/narr]), width=1/20, color=red), k=0..2*narr-1):
q2:=seq( plots:-arrow(eval(M-h*N, [u=0,t=k*Pi/narr]),-eval(d*N, [u=0,t=k*Pi/narr]), width=1/20, color=blue), k=0..2*narr-1):
plots:-display(p1,p2,q1,q2);

My intuition refuses to work about the normals for the Mobius strip.
Mobius strip is not orientable; the normal vector is not continuous. I cannot imagine the "correct" exterior normals for the two "translated" strips.  Isn't the 3D printer going to be confused?

@mmcdara 

It works in Maple 2019. But not properly for Digits>15 (try 25 and 20). Probably Carl will fix this.

Unfortunately you have not included a corresponding Your_NormInv to compare the accuracy.

@mmcdara 

It works in newer versions.
For Maple 2015 just replace z1=z  with  abs(z1-z)<10^(-Digits+1)   (or something similar).

@DJJerome1976 

AFAIK the Risch algorithm is not fully implemented in Maple.
You can see the methods which are tried using
infolevel[int]:=5;

BTW, rewriting just a bit the integrand, Maple integrates easily:

int(sin(x)^(1/3)*(1-sin(x)^2)*cos(x), x);
   

@nm 

Take

ode:=y(x)*diff(y(x),x) - y(x) = 0:
dsolve(ode);

   y(x) = 0, y(x) = x + _C1

This is obviously the correct answer (two solutions; actually there are other solutions by combining these on intervals).

If we solve ode wrt y'  ==>  diff(y(x),x) = 1
[obtained  dividing by y (supposed to be <> 0; that is how solve works)].

Of course, dsolve( diff(y(x),x) = 1 ) ==>  y(x) = x + _C1,
so, y=0 disappears. But y(x)=0 satisfies (of course) the ode.

@Rouben Rostamian  

OP wants an asymptotic expansion for G0(s). He should use asympt(G0(s),s)

The question does not seem to have much sense.

@9009134 

# K[r, s]
for p from 1 to 3 do  
for  s from 1 to 3  do
ff := L[p, s](r, theta, phi)*F[p, 1](r, theta):
Aij:=[indets([F[p,s](r,theta),ff],indexed)[]];
# print(indets=Aij);
g1:=expand(ff, Aij, distributed);
g2:=coeffs(g1, Aij, 'T'):
C1:=int~([g2], theta = 0 .. 2*Pi, r = 0.5 .. 1, epsilon=1e-8, numeric):
kk:=add(C1 .~ T);
kkk := evalf( int(F[p, 1](r, theta)^2, theta = 0 .. 2*Pi, r = 0.5 .. 1) );
k[p, s] := kk/kkk;
print([p, s] = %);
od
 od;

 

@Carl Love 

Unfortunately, as I see, sw_zgeevx_  also fails sometimes.

The bug seems to occur only for multiple (repeated) eigenvalues in the non-symmetric case, so, in practice should be extremely rare.
The examples were found using e.g.

LinearAlgebra:-CompanionMatrix((x-15)^3, x)^+;

 

 

 

@9009134 

1.

f[1,1] is not defined.

2. You must rename my f  to ff  (because you have introduced a f)

ff := L[p, s](r, theta, phi)*F[p, 1](r, theta):
Aij:=[indets(F[p,s](r,theta),indexed)[]];
g1:=expand(ff, Aij, distributed);

3. Now you have B[i,j] etc

Sorry, but I cannot continue with so many changes.


 

@michaelvio 

With series(...,r=r0) you may approximate only near r0.
In your example you have a huge r. Use asympt instead.
[Note that you have not actually used the power of slode; only a few terms, just like the standard series]

restart;
R:=370; Order:=6:
ode:= diff(g(r),r$2)- r/R*g(r)=0:
ic:=g(2*R)=0, D[1](g)(0)=R:
G:=rhs(dsolve([ode,ic],g(r)));

370

 

-740*370^(1/3)*Pi*AiryBi(2*370^(2/3))*AiryAi((1/370)*370^(2/3)*r)/(GAMMA(2/3)*(AiryAi(2*370^(2/3))*3^(2/3)+3^(1/6)*AiryBi(2*370^(2/3))))+740*370^(1/3)*AiryAi(2*370^(2/3))*Pi*AiryBi((1/370)*370^(2/3)*r)/(GAMMA(2/3)*(AiryAi(2*370^(2/3))*3^(2/3)+3^(1/6)*AiryBi(2*370^(2/3))))

(1)

asympt(G, r):

Gapp:=convert(evalf[50](%),polynom);

(0.47817456179569551515317081621201242285756903386552e-602*(1/r)^(1/4)+0.95811159840225010761653791669736639828136520861439e-602*(1/r)^(7/4)+0.14782110748219797338014145566276295493936785779979e-600*(1/r)^(13/4)+0.43638226822222969691007895829878914061306256867382e-599*(1/r)^(19/4))*(exp(19.235384061671344751855362921213373873781867091120/(1/r)^(3/2)))^(1/555)

(2)

evalf[50](eval([G,Gapp], r=400));

[-0.40702431977939449083274686710027843869975497791315e-117, 0.27836877438958079551079281963767137932752590618428e-482]

(3)

evalf[50](eval([G,Gapp], r=1000));

[0.81928332104566902988587362191002268450537068088687e-127, 0.81928332104563623792440214614327162747404927396426e-127]

(4)

evalf[50](eval([G,Gapp], r=10^10));

[0.16056840058550626226900357319617922990909434242673e15051930008898, 0.16056840058550626226900357319617922981745765695502e15051930008898]

(5)

restart;

The integral

J:=Int( exp(I*m*omega + I*b*cos(omega) ), omega=0..2*Pi) ;#          (m integer,  b positive constant)

Int(exp(I*m*omega+I*b*cos(omega)), omega = 0 .. 2*Pi)

(1)

value(J);  # wrong!

0

(2)

 

f:=expand(op(1,J)) assuming m::integer;

exp(I*m*omega)*exp(I*b*cos(omega))

(3)

Prepare the change of variables

g := expand(z^m * exp(I*b*(z+1/z)/2));  z = exp(I*omega);

z^m*exp(((1/2)*I)*b*z)*exp(((1/2)*I)*b/z)

 

z = exp(I*omega)

(4)

G:=simplify(g, [exp(I*omega)=z]) assuming m::integer

z^m*exp(((1/2)*I)*b*z)*exp(((1/2)*I)*b/z)

(5)

answer = 2*Pi*residue(G/z, z=0); # Residue theorem, z=0 essential singularity

answer = 2*Pi*residue(z^m*exp(((1/2)*I)*b*z)*exp(((1/2)*I)*b/z)/z, z = 0)

(6)

fii:=convert(G, Sum, dummy=i);  # we need the constant (coeff of z^0):

z^m*(Sum((((1/2)*I)*b*z)^i/factorial(i), i = 0 .. infinity))*(Sum((((1/2)*I)*b/z)^i/factorial(i), i = 0 .. infinity))

(7)

ai:=eval(op([2,1], fii), z=1); ai_bis:=eval(op([3,1], fii), z=1);

(((1/2)*I)*b)^i/factorial(i)

 

(((1/2)*I)*b)^i/factorial(i)

(8)

Sum(ai * subs(i=i+m, ai), i=0..infinity);  # coeff of z^(-m) in fii

Sum((((1/2)*I)*b)^i*(((1/2)*I)*b)^(i+m)/(factorial(i)*factorial(i+m)), i = 0 .. infinity)

(9)

answer := simplify(2*Pi*value(%)) assuming m::nonnegint;

2*Pi*BesselJ(m, b)*exp(((1/2)*I)*m*Pi)

(10)

evalf(eval([J,answer],[m=3,b=7/3])); #check1

[-0.2869900475e-12-1.168527755*I, -1.168527755*I]

(11)

evalf(eval([J,answer],[m=2,b=5]));   #check2

[-.2925772544-0.3974126502e-12*I, -.2925772544]

(12)

evalf(eval([J,answer],[m=0,b=1]));   #check3

[4.807878861+0.3488081236e-12*I, 4.807878862]

(13)

evalf(eval([J,answer],[m=7,b=4]));   #check4

[-0.2695364704e-12-0.9535405641e-1*I, -0.9535405642e-1*I]

(14)

@Carl Love 

I assumed OP knew what he was doing. The original problem was not stated.

@9009134 

Your constants are not assigned. After giving them numeric values, you may try to compute the integrals, but the integrands being very "odd", you may need to fine-tune int as above.

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