vv

12403 Reputation

19 Badges

8 years, 151 days

MaplePrimes Activity


These are answers submitted by vv

A few considerations.

The equality 
arctan(z/b)/b = int(1/(b^2+u^2), u = 0 .. z)     (1)

is well known when z, b are real and b<>0 (note that I changed the dummy variable, it's not a good idea to use the same letter here).
Your substitution does not make much sense.

Now about the complex case in (1), we must be careful with complex integration. Maple is used for formal (symbolic) computations but it is supposed that the user knows the maths behind these computations.

A complex integral such as (1) is usually a "complex path integral". The integrand is analytic (holomorphic} in C \ {+-I*b}. If the path is missing (as here), it usually means that the integral is independent of the path (i.e. depends only on the endpoints 0 and z). This is true (for our analytic function) if it is considered in a simply connected domain containing 0 and z but not containing the poles +-I*b. So, we must restrict the computations to such domains; for example, C \ {+-I*b} is not simply connected. Ignoring this aspect leads to wrong results or nonsenses.
The conclusion is that you must have knowledge in complex analysis to address such problems.

restart;
Vol:=(u,v,w,U,V,W)->sqrt(-u^2*U^4 + ((-W^2 + u^2 + v^2)*V^2 + (u^2 + w^2)*W^2 - u^4 + (v^2 + w^2)*u^2 - v^2*w^2)*U^2 - v^2*V^4 + ((v^2 + w^2)*W^2 + (v^2 - w^2)*u^2 - v^4 + v^2*w^2)*V^2 - W^2*(W^2*w^2 + (v^2 - w^2)*u^2 - v^2*w^2 + w^4))/12:
m:=20; nr:=0: nrmax:=50;
for U to m do
for V from U to m do
for W from V to U+V-1  do
for u to m do
for v from abs(W-u)+1 to W+u-1  do
for w from max(abs(W-u),abs(U-v))+1 to min(W+u,U+v)-1  do
  vv:=Vol(u,v,w,U,V,W);
  if vv::posint then nr++; lprint(nr, [u,v,w],[U,V,W],vol=vv)  fi;
  if nr>=nrmax then break 6 fi
od:od:od:od:od:od:

# Your example appears at the position 441  (for m=20)

                            m := 20

                          nrmax := 50

1, [8, 8, 7], [2, 3, 4], vol = 6
2, [17, 16, 16], [2, 3, 4], vol = 15
3, [6, 7, 8], [2, 4, 4], vol = 6
4, [6, 8, 7], [2, 4, 4], vol = 6
5, [7, 5, 6], [2, 4, 4], vol = 6
6, [7, 6, 5], [2, 4, 4], vol = 6
7, [4, 8, 7], [2, 5, 6], vol = 6
8, [7, 4, 4], [2, 5, 6], vol = 6
9, [16, 14, 13], [2, 5, 6], vol = 12
10, [4, 6, 5], [2, 7, 8], vol = 6
11, [6, 4, 4], [2, 7, 8], vol = 6
12, [8, 4, 3], [2, 7, 8], vol = 6
13, [12, 14, 13], [2, 7, 8], vol = 24
14, [10, 11, 12], [2, 8, 8], vol = 21
15, [10, 12, 11], [2, 8, 8], vol = 21
16, [11, 9, 10], [2, 8, 8], vol = 21
17, [11, 10, 9], [2, 8, 8], vol = 21
18, [17, 11, 12], [2, 8, 8], vol = 21
19, [17, 12, 11], [2, 8, 8], vol = 21
20, [8, 12, 11], [2, 9, 10], vol = 21
21, [11, 8, 8], [2, 9, 10], vol = 21
22, [8, 10, 9], [2, 11, 12], vol = 21
23, [10, 8, 8], [2, 11, 12], vol = 21
24, [17, 8, 8], [2, 11, 12], vol = 21
25, [12, 8, 7], [2, 13, 14], vol = 24
26, [16, 6, 5], [2, 13, 14], vol = 12
27, [16, 22, 21], [2, 13, 14], vol = 60
28, [17, 3, 4], [2, 16, 16], vol = 15
29, [17, 4, 3], [2, 16, 16], vol = 15
30, [12, 19, 19], [2, 19, 19], vol = 72
31, [8, 24, 23], [2, 19, 20], vol = 42
32, [17, 16, 17], [3, 5, 7], vol = 30
33, [4, 7, 8], [3, 6, 6], vol = 9
34, [4, 8, 7], [3, 6, 6], vol = 9
35, [14, 16, 17], [3, 6, 6], vol = 36
36, [14, 17, 16], [3, 6, 6], vol = 36
37, [4, 6, 6], [3, 7, 8], vol = 9
38, [7, 11, 10], [3, 7, 8], vol = 24
39, [8, 4, 2], [3, 7, 8], vol = 6
40, [11, 7, 8], [3, 7, 8], vol = 24
41, [16, 16, 17], [3, 7, 8], vol = 48
42, [16, 8, 10], [3, 8, 10], vol = 21
43, [11, 15, 16], [3, 9, 9], vol = 42
44, [11, 16, 15], [3, 9, 9], vol = 42
45, [7, 8, 7], [3, 10, 11], vol = 24
46, [9, 15, 14], [3, 11, 12], vol = 48
47, [8, 14, 13], [3, 13, 14], vol = 48
48, [9, 12, 11], [3, 14, 15], vol = 48
49, [9, 17, 15], [3, 15, 16], vol = 42
50, [11, 9, 9], [3, 15, 16], vol = 42

Most "complex plot" commands are very simple, but are convenient. Basically they just take the real and imaginary parts and call a regular plot.
For example, for f:=sin(x+I):

complexplot(f,x=-Pi..Pi);

is equivqlent to:

plot([Re(f),Im(f),x=-Pi..Pi]);
restart;
i2 := (x,y) -> -(1/2)*I*(exp(I*x)*(sin(x)/x)-exp(I*y)*(sin(y)/y))/(x-y);
i3_r := -(1/2)*I*(i2(y,z)-i2(y,x))/(z-x);
normal(eval(series(eval(i3_r,[x=t*x,y=t*y,z=t*z]),t=0),t=1));

On my rather old computer, for any of your graphs (having 21 vertices and 76 edges), VertexConnectivity needs about 0.5 sec.
So, it is easy to estimate the total time for your almost 15000 graphs.

It should be easy to write a compiled version for VertexConnectivity. It will be much faster (and a skilled user can do it), but probably the designers did not anticipate that someone will need it for so many graphs.

`simplify/size` is responsible for this. It calls `simplify/size/size` to compute the "math complexity" (see ?simplify/size) and chooses the simpler one. 

ex:=[t^2*x^2+t^2*y^2, t^2*x^2-t^2*y^2]:
`simplify/size/size`~(%);
#         [56, 59]

simplify~(ex);
#         [t^2*(x^2+y^2), t^2*x^2-t^2*y^2]

`simplify/size/size`~(%);
#        [38, 59]

Note that simplify(..., size=false) keeps both expressions unchanged.

If you want the rational exponent of x^r, where x is a variable and r is rational (including r=1)
then r=1 must be treated separately.

Exponent := proc(x::{name, name ^ rational}) 
  if x::name then 1 else op(2, x) fi 
end proc:
Exponent(x^(-3/2));   # -3/2
Exponent(Y[1,3]);     # 1

Note that r=0 is not possible because x^0 is automatically simplified to 1.

remove removes operands (1st level) of an expression.
x has one operand (x itself)
diff(y(x), x) has two operands (y(x) and x)  [so, you can remove only y(x) or x or both].
2*diff(y(x), x) has two operands (2 and diff(y(x), x))

Wind:=proc(L::listlist, P::list, isclosed:=true)
  local p,u,k, x0:=P[1],y0:=P[2];
  if member(P,L) then return undefined fi;
  u:=seq(arctan(p[2]-y0,p[1]-x0),p=L),  `if`(isclosed=true,arctan(L[1][2]-y0,L[1][1]-x0), NULL);
  (u[-1]-u[1])/(2*Pi) + add(piecewise( u[k+1]-u[k]<-Pi,+1, u[k+1]-u[k]>Pi,-1, abs(u[k+1]-u[k])=Pi, undefined, 0), k=1..nops([u])-1)
end:

It returns 0 if the point is outside, and undefined if it is on the boundary.

L:=[[0,0],[200,0],[200,300],[250,400], [500,300],[500,500],[0,500]]:
Wind(L, [300,300]), Wind(L, [100,200]), Wind(L, [0,50]);

                        0, 1, undefined

N depends on the order in which the numbers are eliminated.
The largest value of N is N_max = 2^(n-1) * (3*n - 4) + 2 
(you have n = 2012).
Try your hand for the smallest value!

It's obvious that the smartview option cannot guess all user's intentions.
If you don't know the relevant y-domain (used by view=...) you may try a slope criterion. E.g.

restart;
f:=exp(-3^(1/2)*(cos(x)-1)/sin(x)):
f1:= diff(f, x):
maxslope:=100:
ff:=piecewise(abs(f1) < maxslope, f, undefined):
plot(ff, x=0..2*Pi);

Maple can compute the integral directly as a hypergeom function.
So, this is for didactic purposes only.

 

 

restart;

with(IntegrationTools):

J(n)=Int(1/(x^2+1)^n, x):

Parts(%, 1/(x^2+1)^n):

eval(%, x^2=q(x)-1):

simplify(Expand(%)) assuming n::integer:

eval(%, [ Int(q(x)^(-n), x)=J(n), Int(q(x)^(-n-1), x)=J(n+1) ]):

solve(%, J(n+1)): J(n+1) = simplify(eval(%, q(x)=x^2+1)); # the recurrence relation

J(n+1) = (1/2)*(x*(x^2+1)^(-n)+J(n)*(2*n-1))/n

(1)

R:=unapply(rsolve({%, J(1)=arctan(x)}, J(n)), n);

proc (n) options operator, arrow; (1/2)*GAMMA(n-1/2)*(Sum(x*Pi^(1/2)*GAMMA(n1+1)/(n1*(x^2+1)^n1*GAMMA(n1+1/2)), n1 = 1 .. n-1)+2*arctan(x))/(Pi^(1/2)*GAMMA(n)) end proc

(2)

value(R(6));

(63/256)*x/(x^2+1)+(21/128)*x/(x^2+1)^2+(21/160)*x/(x^2+1)^3+(9/80)*x/(x^2+1)^4+(1/10)*x/(x^2+1)^5+(63/256)*arctan(x)

(3)

simplify(diff(%,x)); #check

1/(x^2+1)^6

(4)

 

restart;
Digits := 15: 
nu_0 := sqrt(k^2 - k_0^2): nu_m := sqrt(k^2 - k_m^2): 
pref2 := 2*rho_me*exp(-nu_0*(zs + z))*BesselJ(0, k*sqrt((xs - x)^2 + (ys - y)^2))*k/(rho_me*nu_0 + rho_0*nu_m*tanh(nu_m*d)):

f:= (eval(0.5*(((eval(pref2, [z = 0, k_0 = 0.018371886863098, xs = 0, ys = 0, zs = 10, rho_0 = 1.2, d = 0.04, rho_me = 1.528516816439260 - 1235.297048886680*I, k_m = 0.490806242885258 - 0.490314205803914*I])))), [z = 0, k_0 = 0.018371886863098, xs = 0, ys = 0, zs = 10, rho_0 = 1.2, d = 0.04, rho_me = 1.528516816439260 - 1235.297048886680*I, k_m = 0.490806242885258 - 0.490314205803914*I]) ):

dom := [k=0..100, x = -2/2 .. 2/2, y = -3/2 .. 3/2]:  # 100 is enough
evalf(Int(f, dom, digits=10));

                               0.5867621705 - 0.1097638086*I

The answer is correct, a being a name and non-constant.
eval(diff(y(x,t), x$3), x=Pi);  and  eval(diff(y(x,t), x$3), x=a+1);

will look as you expect.

You may declare a as a constant. In this case the result will appear as for Pi:
constants := constants, a;
eval(diff(y(x,t), x$3), x=a)
;
            

restart;
N:=100:  dx:=evalf(2*Pi/N):
X:=<seq(k*dx, k=-N/2..N/2)>:
Y:= sin~(X):
f:=unapply(CurveFitting:-Spline(X, Y, x, degree=1),x):
Z:= [seq( evalf(Int(f, X[1]..x, epsilon=1e-5)), x=X)]:
plots:-display(plot(X,Y),  plot(X,Z), color=[red,blue]); 

First 11 12 13 14 15 16 17 Last Page 13 of 111