## dharr

Dr. David Harrington

## 6337 Reputation

19 years, 352 days
University of Victoria
Professor or university staff

## Social Networks and Content at Maplesoft.com

I am a professor of chemistry at the University of Victoria, BC, Canada, where my research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

## arrows palette...

The arrows palette has diagonal arrows, which can be inserted in 2D, e.g. in a Table:

 -5 -2 0 2 3 4 | - 0 + 0 - 0 - 0 + | | a lok. min b lok. max a van v.tg a lok. min b |

Here's one way.

 >
 >
 >
 >

 >

## Mathcontainer...

The font size in a TextArea component cannot be changed as far as I can see. But if you change to a MathContainer component, you can use all the typesetting capabilities. Here is one way to do it:

S4ThalèsAnimation.mw

## Zeckendorf's theorem...

Zeckendorf's theorem in the title of your post is about a unique sum of Fibonacci's, not all sums. @mmcdara has just given a version of all sums; here is an efficient implementation of the unique version.

Zeckendorf's theorem. There is a unique sequence of distinct fibonacci numbers adding to F if we exclude consecutive fibonacci numbers.

 >

Golden ratio

 >

n for largest fibonacci not greater than F - see wikipedia

 >

Largest fibonacci not greater than 100.

 >

Greedy algorithm - successively find largest fibonacci's f not greater than the remainder r.

 >

 > zeck := proc(n::posint)   local nlarge, r, f, i;   nlarge := F -> floor(log[(1 + sqrt(5))/2](sqrt(5)*(F + 1/2)));   r := n;   f := table();   for i do     f[i] := combinat:-fibonacci(nlarge(r));     r := r - f[i];   until r = 0;   convert(f, list); end proc:
 >

 >

 >

## zero outside range...

Using trace and showstat shows that this doesn't go through GAMMA, despite the help page. Hope I've traced this correctly; I haven't given all the details. (For positive integers, binomial(a,b) = 0 for b>a is consistent with no ways to choose b out of a.)

 > restart;
 > trace(binomial);

Consider first integers outside the range 0<=k<=n; the result is just zero

 > binomial(4, 6);

execute binomial:-ModuleApply, args = 4, 6

{--> enter binomial:-ModuleApply, args = 4, 6

<-- exit binomial:-ModuleApply (now at top level) = 0}

The code is just:
elif type(A,'integer') and type(B,'integer') then
if 0 < A then
if B < 0 or A < B then
res := 0 (result)

Now consider

 > binomial(n, n + 2) assuming n::posint;

execute binomial:-ModuleApply, args = \`n~\`, n+2

{--> enter binomial:-ModuleApply, args = \`n~\`, n+2

<-- exit binomial:-ModuleApply (now in \`assuming\`) = 0}

So what happens here:

B:=normal(b) -> n+2;
Work out AmB: = a-B -> -2
A:=normal(a) -> n;
s: = signum(A)  -> 1

If s = 1 or -1 and AmB::negint (yes)

then

i := signum(0,B,0) -> 1 (sign of B treating 0 as 0)
if i is 1 or -1 then
if s = 1 or s = i then
res := 0

So result is zero.

In other words, since A and B are positive and A-B is a negative integer the result is zero

binomial(n, n+1) uses  `tool/type` which makes it a bit more obscure but seems to be a similar argument.

 >

This is not as systematic as @mmcdara's method, but answers the question for your second example containing sigma__d3. You could use the approximate form for Lambda instead.

 >
 > local gamma:with(plots):

# Define the assumptions - or Maple wouldn't know

 > assume(0 < gamma, 0 < sigma__d, 0 < sigma__d3, 0 < sigma__v); interface(showassumed=0);

# Input lambda parameters

 > Lambda := RootOf(8*_Z^4 + 12*Gamma*_Z^3 + (5*Gamma^2 - 4)*_Z^2 - 4*Gamma*_Z - Gamma^2); f__1 := Lambda*(sigma__v/sigma__d): f__2 := f__1: f__3 := sqrt(2)*sigma__v/sigma__d3:

Let's scale these. If we divide f__1 by sqrt(2)*(sigma__v/sigma__d) we get Lambda/sqrt(2); we divide f__2 and f__3 by the same.

 > f__1sc := f__1/(sqrt(2)*(sigma__v/sigma__d)); f__2sc := f__2/(sqrt(2)*(sigma__v/sigma__d)); f__3sc := f__3/(sqrt(2)*(sigma__v/sigma__d));

Now Gamma = gamma*sigma__v*sigma__d = gamma*sigma__v*sigma__d3*f3sc = a*f3sc, where a = gamma*sigma__v*sigma__d3.

We want to know when f3sc = f1sc, or Gamma/a = Lambda/sqrt(2), which will only depend on the combined parameter a = gamma* sigma__v*sigma__d3

 > Gamma/a=f__1sc; aval:=solve(%,a);

Plot f__1sc = f__3sc = sigma__d/sigma__d3 vs a = gamma*sigma__v*sigma__d3

 > pp:=plot([aval,f__1sc,Gamma=0..2],labels=[gamma*sigma__v*sigma__d3,sigma__d/sigma__d3]):
 > display(pp,textplot([3,0.3,typeset('f__3' < 'f__1')]),textplot([3,0.45,typeset('f__3' > 'f__1')]));

Test for some parameters

 > params:={gamma = 1, sigma__v = 1, sigma__d3 =3, sigma__d = 1}; evalf(eval([gamma*sigma__v*sigma__d3, sigma__d/sigma__d3], params));

This point is below the line, so we expect f3 < f1

 > evalf(eval(f__1sc, Gamma=eval(gamma*sigma__v*sigma__d,params))); evalf(eval(f__3sc, params));

 >

## results can be different...

The documentation for IsFrobeniusGroup gives an example where IsFrobeniusGroup is true but IsFrobeniusPermGroup is false, so I think this result is OK. The documentation you referred to (for FrobeniusGroup) refers to the situation where  IsFrobeniusPermGroup is true, and so IsFrobeniusGroup will be true.

## different perm sets...

The small group (8,3) is expressed as a permutation group on 1,2,..8, but DihedralGroup(4) is expressed as a permutation group on 1,2,3,4. So the domain [1,2,3,4] is only half the points for (8,3) but it is all the points for D4. ReducedDegreePermGroup will convert (8,3) to the form you need.

 >
 >
 >

 >

 >

 >

 >

 >

 >

 >

 >

 >

 >

## congruent numbers...

You have missed the fact that the algorithm is only for squarefree integers. For example 20 should give true, but to decide this you need to divide by the largest square factor 4 and apply the algorithm to the quotient 5. You also missed 6, which is squarefree, for some reason that wasn't clear to me.

(I used @Ronan's version replacing print("True") with return true etc, but have my own working version if you are interested.)

## numerical issues...

As @C_R and @Rouben Rostamian have pointed out, the function is not specified accurately enough.

 > restart;
 > with(plots):
 > F:=phi->3.924999 - 0.024999/sqrt(1 - 2*phi) - 3.900/(1 - phi/6)^(3/2) - 1.648094618*10^(-14)*sqrt(3)*sqrt(1836)*(((1.3972 + sqrt(3)*sqrt(1836))^2 + 3672*phi)^(3/2) - (1.3972 + sqrt(3)*sqrt(1836))^3 - ((1.3972 - sqrt(3)*sqrt(1836))^2 + 3672*phi)^(3/2) + ((1.3972 - sqrt(3)*sqrt(1836))^2)^(3/2))-(1/18)*(sqrt(3)*sqrt(300)*(((1.4472 + sqrt(3)*sqrt(300)/300)^2 - 2*phi)^(3/2) - (1.4472 + sqrt(3)*sqrt(300)/300)^3 - ((1.4472 - sqrt(3)*sqrt(300)/300)^2 - 2*phi)^(3/2) + (1.4472 - sqrt(3)*sqrt(300)/300)^3)):

Problem seems numerically poorly posed: value at zero  (intended to be zero?) depends on floatng point 3.924999 - 0.024999 = 3.900, and then has term with coefficient 1.648e-14; large values also (3672)

Function goes to (close to?) zero at 0 and near 0.095

 > plot(F(phi),phi=-0.2..0.2,-0.00001..0.00001);

So the integrand goes to infinity/large value at these points.

 > integrand:=1/sqrt(-2*F(phi)):

Unreliable evaluation at default digits (see Roubens plot)

 > evalf(eval(integrand,phi=1e-9)); evalf(eval(integrand,phi=1e-8));

 > plot(integrand,phi=-1..1);

 > plot(integrand,phi=-1.8..-0.4);

Complex

 > evalf(eval(integrand,phi=-1.6));

Rather than repeatedly solving integrals in a plot, can pose it as an ode, Arbitrarily integrate from -1.2

 > ode:=diff(x(phi),phi)=integrand:
 > ans1:=dsolve({ode,x(-1.2)=0.},numeric);

Fails close to zero while still increasing, suggestng that the area might be infinite

 > p1:=plots:-odeplot(ans1,[phi,x(phi)],-1.2..0.5): display(p1);

Warning, cannot evaluate the solution further right of -.97056569e-3, probably a singularity

Restart integration on the other side of zero

 > ans2:=dsolve({ode,x(1e-2)=0},numeric);

 > p2:=plots:-odeplot(ans2,[phi,x(phi)],1e-2..0.4): display(p2);

 >

## TravelingSalesman for digraph...

Although you have to return to the starting point, you can force it to return via the last vertex by using a directed graph, and so find the shortest path between the first and last vertices. I'm not sure I fully understood your description, so this may not be exactly what you want. I worked out edge distances from coordinates, but you could just enter edge weights for each edge.

 >
 >

We want shortest distance between first and last vertex (1 and 5 here) that visits all vertices
List of coordinates of points labeled 1..n

 >

Edges between points - include edge {1,5} between start and end vertices

 >

Find geometric distances along edges

 >

Ensure return is via vertex 5 and not any other neighbors of vertex 1

 >
 >

 >

 >

 >

 >

 >

## without a literature result...

Here's a way that avoids using the known result for the othogonality. On the other hand it should be easier than finding a generating function...

 > restart; with(PDEtools): #with(Maplets[Elements]):
 > eq0:= E1 = (n+1/2)*_h; param := m1 =_h*sigma^2;

 > os := expand(isolate((-_h^2/(2*m1))*diff(psi(x),x,x)+(m1*x^2/2)*psi(x)-E1*psi(x)=0,diff(psi(x),x,x)));

 > eq1 := PDEtools:-dchange([x=sqrt(_h/m1)*zeta,psi(x)=f(zeta)*exp(-zeta^2/2)],os,[zeta,f(zeta)],params={m1,_h,E1}); eq2 := expand(isolate(eq1,diff(f(zeta),zeta\$2))); eq3 := simplify(eval(remove(has,convert(rhs(dsolve(eq2,f(zeta))),Hermite),KummerM),eq0)*exp(-zeta^2/2),radical)       assuming zeta::positive; eq4 := simplify(subs(param,simplify(subs(zeta=sqrt(_h/m1)*x,eq3)))) assuming sigma::positive; eq5 := unapply(eq4,n,x,sigma);

Work out the first few integrals

 > ints:=[seq(int(simplify(eq5(n,x,sigma)^2,'HermiteH'),x=-infinity..infinity),n=0..10)] assuming positive;

Find out the pattern

 > gfun:-guessgf(ints,x); gf:=%[1]: series(gf,x,11): [op(convert(%,polynom))]: seq(%[i]*(i-1)!,i=1..nops(%));

The nth coefficient of the exponential generating function is the nth derivative evaluated at x=0
So the general integral is

 > eval(diff(gf,[x\$n]),x=0); intn:=convert(simplify(%),factorial) assuming n::nonnegint;

And so the coefficient is given by

 > eqc2:=c__2=solve(intn=1,c__2)[1];

Leading to the final result for the wavefunction

 > simplify(eval(eq4,eqc2)) assuming n::nonnegint: unapply(%,n,x,sigma);

 >

## 0.61822293776...

I can only get 11 place accuracy without going beyond hardware precision. Probably a better estimate of the error term is required.

 > restart;

We want the integral of the absolute value of the following from 0 to infinity to 14 digits of precision

 > y:=sin(x^4)/(sqrt(x) + x^2);

 > plot(abs(y),x=0..3);

We could integrate the successive "bumps" and sum them out to a large number of bumps, but convergence is slow and the error of the omitted area is hard to estimate.
Consider a bump over the range X[i] to X[i+1] for some large integer i

 > X := i->(Pi*i)^(1/4);

 > dx:=X(i+1)-X(i); midx:=(X(i)+X(i+1))/2; I1:=Int(y,x=X(i)..X(i+1));

The bump approximates a scaled sin function for large i, with height 1/(sqrt(x) + x^2) and width dx. Even for i=10 this is visually a good approximation

 > plot(eval([[y,1/(sqrt(midx)+midx^2)*sin((x-X(i))/dx*Pi)],x=X(i)..X(i+1)][],i=10),color=[red,blue],tickmarks=[default\$2]);

The area of a half period of a sine function is 2, and the area of the rectangle enclosing it is Pi, so the area of a bump is approximately

 > A:=2/Pi*dx/(sqrt(midx) + midx^2):

Compare with "exact" value for some i values - show absolute and relative errors

 > evalf[20](eval([Int(y,x=X(i)..X(i+1)),A],i=5000)); %[2]-%[1],(%[2]-%[1])/%[1]; evalf[20](eval([Int(y,x=X(i)..X(i+1)),A],i=6000)); %[2]-%[1],(%[2]-%[1])/%[1];

So the area from some large i to infinity is approximately 2/Pi times the integral of the envelope function.
The relative error in this should be better than the relative error above for the same i, so for 6000 the relative error is about 1e-9, or correct to about 10 decimal places

 > Int(2/Pi/(sqrt(x) + x^2),x=X(i)..infinity); evalf[20](eval(%,i=6000)); %*1e-9;

So numerically integrate the first n bumps and sum (backwards to reduce loss of significance), then add the above estimate for the rest.
At Digits=15 (hardware precision) we have n=1000 gives , n=5000 gives , n=6000 gives

At Digits:=17 (slow), n=1000 gives , suggesting the Digits=15 numerical integrations are in error only in the last digit, i.e. we have 14 digits of prcesion.
However estimate from n=5000 to n=6000 still changing in the 12th place due to the error term, slightly better than as expected from above

 > Digits:=15; n:=6000; # n even seq[reduce=`+`](evalf(eval(I1,i=ii)),ii=n-2..0,-2) +seq[reduce=`+`](evalf(eval(-I1,i=ii)),ii=n-1..1,-2) +evalf(Int(2/Pi/(sqrt(x) + x^2),x=X(n)..infinity));

Going to more, say 10000, terms will degrade the numerical sum and not much improve the error term, but probably the best possible with this strategy and hardware precision.
n = 10000 gives , n=20000 gives

 > Digits:=15;#Digits:=17; n:=20000; # n even seq[reduce=`+`](evalf(eval(I1,i=ii)),ii=n-2..0,-2) +seq[reduce=`+`](evalf(eval(-I1,i=ii)),ii=n-1..1,-2) +evalf(Int(2/Pi/(sqrt(x) + x^2),x=X(n)..infinity));

 >

## Same results...

When you say they are not consistent, are you suggesting they are not mathematically equivalent? If so do you have some reason to think that? The following indicates they are the same, though it is frustrating that simplify doesn't show that directly. (Following is Maple 2023.)

 >
 >
 >
 >

 >

This should give zero but it doesn't.

 >

This works but it assumes all indeterminates are real, which is not what you want. Nonetheless it suggests they are probably equal.

 >

Choose some random complex values for the indeterminates and check

 >

 >

## Maple can tell you...

Of course the more specific you are about the function h(x) and the initial/boundary condition, the more detailed the solution will be. The help page is found online here.

 >
 >

 >

 >

Following gives a help page describing the Bernoulli case and the transformation used to solve it.

 >