## Not all the zeros of LegendreP are computed...

Hi,

Working with Legendre Polynomials (LegendreP) I observed that solve doesn't find the correct number of zeros.
More precisely, for N > 17, solve(LegendreP(N, x)) finds less zeros than N.

I wrote a procedure based on a theorem about the intertwined location of the zeros of orthogonal polynomial of successive degrees. So this problem is not blocking, but I would like to understand while solve(LegendreP(N, x)) doesn't always do the job.

 > restart:
 > Z := n -> op~(2, { allvalues(solve(LegendreP(n,x))) } );
 (1)
 > Digits:=10: Z(17): numelems(%);
 (2)
 > Z(18): numelems(%);
 (3)
 > Digits:=15: Z(18): numelems(%);
 (4)
 > Digits:=20: Z(18): numelems(%);
 (5)
 > Zf := n -> op~(2, { allvalues(solve(evalf(LegendreP(n,x)))) } ); Z(18): numelems(%);
 (6)
 > # Let z[N][i] the ith zero of any orthogonal polynomial P(N,x) of degree N. # # It is known that each open interval(z[N][i], z[N][i+1]) contains # exactly a unique zero of the of P(N+1,x). Z17 := [ -1, Z(17)[], 1]: Z18 := NULL: for n from 1 to 18 do   Z18 := Z18, fsolve(LegendreP(18,x),  x=Z17[n]..Z17[n+1]); end do: numelems({Z18})
 (7)
 > # A procedure to compute zeros of LegendreP up to degree N zeros := proc(N)   local zeros_table, Z, n, p, z:   zeros_table := table([0=[]]):   Z := [-1, 1]:   for n from 1 to N do     z := NULL:     for p from 1 to n do       z := z, fsolve(LegendreP(n,x),  x=Z[p]..Z[p+1]);     end do;     zeros_table[n] := [z]:     Z := [-1, z, 1]   end do;   return zeros_table end proc:

## Why dsolve (numeric, maple2015) stop calculating ...

I used dsolve to solve the Initial value problems numericaly.
When I set the parameter range=1..5*10^4 , it works and cost only about 200s cpu time.
But if I set range=1..2*10^5, it stop running ( cpu time stop) when the mserver memory reach about 1.5 G. (the memory record in the bottom-right of maple interface is about 700M .)

## Why is precedence losts and how to recover it?...

Hi,

Using InertfForm, I found this strange result:

 > e := (A+B)*(A+C); ife  := InertForm:-Parse(convert(e, string)); InertForm:-Display(ife);
 (1)
 >

Why doesn't Display(..) return (A+B)*(A+C) ?

## How to manage the number of colors colorscheme use...

Hi,

I often need to draw several curves on the same figure and the question of which colors to choose to distinguish them is crucial.
In particular when these curves depend on a parameter I would like to color them by using a "regular" change of color.
The trick I use is this one

```# here I want to plot N=10 different curves

N := 10:

# step 1: generate a list of colours using colorscheme
# (I was hoping for 10 levels of colors but I got 19)

plot(x, x=0..1, colorscheme=["Gold", "Blue"], numpoints=N):
MyColors := op([1,2, 2], %);
M := numelems(MyColors[..,1]);

# step 2: do the plots of interest, for instance

plots:-display(
seq(
plot(m*x^2, x=0..1, color=ColorTools:-Color([entries(MyColors[round(m/N*M)], nolist)])),
m=1..N
)
);
```

My question is: How can we make that MyColors be a matrix whose the number of rows equal N (this to avoid the round(m/N*M) operation)?
(I will accept any other strategy)

## A digression about the US electoral college tie

by: Maple 2015

A fascinating race is presently running (even if the latest results seem  to have put an end to it).
I'm talking of course about the US presidential elections.

My purpose is not to do politics but to discuss of a point of detail that really left me puzzled: the possibility of an electoral college tie.
I guess that this possibility seems as an aberration for a lot of people living in democratic countries. Just because almost everywhere at World electoral colleges contain an odd number of members to avoid such a situation!

So strange a situation that I did a few things to pass the time (of course with the earphones on the head so I don't miss a thing).
This is done with Maple 2015 and I believe that the amazing Iterator package (that I can't use thanks to the teleworking :-( ) could be used to do much more interesting things.

 > restart:
 > with(Statistics):
 > ElectoralCollege := Matrix(51, 2, [
 > Alabama,        9,        Kentucky,        8,        North_Dakota,        3,
 > Alaska,        3,        Louisiana,        8,        Ohio,        18,
 > Arizona,        11,        Maine,        4,        Oklahoma,        7,
 > Arkansas,        6,        Maryland,        10,        Oregon,        7,
 > California,        55,        Massachusetts,        11,        Pennsylvania,        20,
 > Colorado,        9,        Michigan,        16,        Rhode_Island,        4,
 > Connecticut,        7,        Minnesota,        10,        South_Carolina,        9,
 > Delaware,        3,        Mississippi,        6,        South_Dakota,        3,
 > District_of_Columbia,        3,        Missouri,        10,        Tennessee,        11,
 > Florida,        29,        Montana,        3,        Texas,        38,
 > Georgia,        16,        Nebraska,        5,        Utah,        6,
 > Hawaii,        4,        Nevada,        6,        Vermont,        3,
 > Idaho,        4,        New_Hampshire,        4,        Virginia,        13,
 > Illinois,        20,        New_Jersey,        14,        Washington,        12,
 > Indiana,        11,        New_Mexico,        5,        West_Virginia,        5,
 > Iowa,        6,        New_York,        29,        Wisconsin,        10,
 > Kansas,        6,        North_Carolina,        15,        Wyoming,        3 ]):
 (1)
 (2)
 > ec := convert(ElectoralCollege, listlist):
 > # Sets of states that form an electoral college tie R      := 10^5: nbties := 0: states := NULL: for r from 1 to R do   poll  := combinat:-randperm(ec):   cpoll := CumulativeSum(op~(2, poll)):   if tie in cpoll then     nbties := nbties+1;     place  := ListTools:-Search(tie, cpoll);     states := states, op~(1, poll)[1..place]:   # see below   end if: end do:
 > # electoral college tie is not so rare an event # (prob of occurrence about 9.4 %). # # Why the hell the US constitution did not decide to have an odd # number or electors to avoid ths kind of situation instead of # introducing a complex mechanism when tie appears???? nbties; evalf(nbties/R); states := [states]:
 (3)
 > # What states participate to the tie? names := sort(ElectoralCollege[..,1]): all_states_in_ties := [op(op~(states))]: howoften := Vector(                     51,                     i -> ListTools:-Occurrences(names[i], all_states_in_ties)             ): ScatterPlot(Vector(51, i->i), howoften);
 > # All the states seem to appear equally likely in an electoral college tie. # Why? Does someone have a guess? # # The reason is obvious, as each state must appear in the basket of a candidate, # then in case of a tie each state is either in op~(1, poll)[1..place] (candidate 1) # or either in op~(1, poll)[place+1..51] (candidate 2); # So, as we obtained 9397 ties, each states appears exactly 9397 times (with # different occurences in the baskets of candidate 1 and 2).
 > # Lengths of the configurations that lead to a tie. # # Pleas refer to the answer above to understand why Histogram(lengths) should be # symmetric. lengths := map(i -> numelems(states[i]), [\$1..nbties]): sort(Tally(lengths))
 (4)
 > Histogram(lengths, range=min(lengths)..max(lengths), discrete=true)
 > ShortestConfigurations := map(i -> if lengths[i]=min(lengths) then states[i] end if, [\$1..nbties]): print~(ShortestConfigurations):
 (5)
 > LargestConfigurations := map(i -> if lengths[i]=max(lengths) then states[i] end if, [\$1..nbties]): print~(LargestConfigurations):
 (6)
 > # What could be the largest composition of a basket in case of a tie? # (shortest composition is the complementary of the largest one) ecs   := sort(ec, key=(x-> x[2])); csecs := CumulativeSum(op~(2, ecs)): # Where would the break locate? tieloc := ListTools:-BinaryPlace(csecs, tie); csecs[tieloc..tieloc+1]
 (7)
 > # This 40  states coniguration is not a tie. # # But list all the states in basket of candidate 1 and look to the 41th state (which is # in the basket of candidate 2) ecs[1..tieloc]; print(): ecs[tieloc+1]
 (8)
 > # It appears that exchanging Virginia and New_Jersey increases by 1 unit the college of candidate 1 # and produces a tie. LargestBasketEver := [ ecs[1..tieloc-1][], ecs[tieloc+1] ]; add(op~(2, LargestBasketEver))
 (9)
 > # The largest electoral college tie contains 40 states (the shortest 11)

## How to plot the difference between 2 dsolve[numer...

Hi,

I have an ODE system with a free parameter P, that I can only solve numerically.
To simplify let's suppose this parameter has only 2 values P1 and P2.
Let SOL1 the solution of this system when P=P1 and SOL2 its solution for P=P2:

Again for a sake of simplification let's suppose the independent variable is t and the dependent one is x(t).

How can I plot the difference x1(t)-x2(t) where x1(t) is the solution corresponding to SOL1 and x2(t) the solution corresponding to SOL2
(sorry I screwed something up and changed the config of my keyboard making the interrogation point unaccessible)

Here is a detailed example

 > restart:
 > with(plots):
 > # This is a 1D problem. # A body with unit mass is submitted to an acceleration and to a solid-solid friction resistance # The friction model is a simple Coulomb's sys := {          diff(x(t), t) = v(t),          diff(v(t), t) = t - piecewise(v(t) > 0, f, v(t) < 0, -f, 0),          x(0) = 0,          v(0) = 0        }; sol__0 := dsolve(sys, numeric, parameters=[f], maxfun=10^6)
 (1)
 > # This system doesn't seem to be solved formally sol__exact := dsolve(sys, {v(t), x(t)})
 > # Attempt of solving it numerically. # # Let's observe that the "v" solution must be an increasing function of t # so the friction model could be replaced by piecewise(v(t) > 0, f, 0) for v(t) >=0. # But this is a particular case (think of a sine acceleration for instance) and # I keep using the initial Coulomb's model. # # The main problem, as it appears below is the discontinuity of the friction effect at v(t)=0 # (thus at the initial time). sol := dsolve(sys, numeric, parameters=[f], maxfun=10^6): sol(parameters=[1]): sol(0.1);
 > # Clearly increasing maxfum will help to simulate on a longer time but it's not the # most efficient way to handle this problem. # # The classical way is to use a regularized Coulomb's model ; for instance sys := {              diff(x(t), t) = v(t),              diff(v(t), t) = t - f *tanh(v(t)/c),              x(0) = 0,              v(0) = 0            }; sol__1 := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8)
 (2)
 > # The smaller c the closer the solution to the one invoking the Coulomb's model. # # The question is: how can we measure as close to sol__exact we are? # The idea is to choose a "small" value of c, lets say c__0, and to declare that the solution # obtained with c__0 is "the exact solution". # # Studying the convergence of the solution as c__0 becomes smaller and smaller  will help # choosing c__0. # Of course the computational time is an increasing function of 1/c__0 : the smaller c__0 # the higher this time. c__0 := 1e-8:  # arbitrary choice sol__ref := dsolve(sys, numeric, parameters=[f, c], maxfun=10^8): sol__ref(parameters=[1, c__0]): tmax := 10: CodeTools:-Usage(sol__ref(tmax)): numfun__0 := sol__ref(numfun); odeplot(sol__ref, [t, v(t)], t=0..tmax)
 memory used=33.20KiB, alloc change=0 bytes, cpu time=34.76s, real time=36.19s, gc time=0ns
 > # Assuming c__0 = 1e-8 is a value small enough to consider the solution plotted above # is "the solution for the true Coulomb's model". # # Here are a few solutions for larger values of c: cs   := 10^~[\$-6..-1]: tmax := 10: printf("c = 1e-08  number of calls = %d\n", numfun__0): for i from 1 to nops(cs) do   sol__1(parameters=[1,cs[i]]):   sol__1(tmax):   couleur := ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]):   Xplot__||i := odeplot(                          sol__1, [t, x(t)], t=0..tmax,                          legend=cs[i],                          color=couleur                        ):   printf("c = %1.0e  number of calls = %d\n", cs[i], sol__1(numfun)) end do: LY  := plottools:-transform((x, y) -> [x, log(y)]): LXY := plottools:-transform((x, y) -> [log(x), log(y)]): display(          seq(Xplot__||i , i=1..nops(cs)),          view=[0..3, 0..2], labels=[t, x], title="x(t) (lin-lin)"        );
 c = 1e-08  number of calls = 50831166 c = 1e-06  number of calls = 1308504 c = 1e-05  number of calls = 126848 c = 1e-04  number of calls = 12825 c = 1e-03  number of calls = 1513 c = 1e-02  number of calls = 456 c = 1e-01  number of calls = 246
 > # From this plot it seems that values of c less than 1e-4 give roughly the same curves. # # Remembering that the cimputational time is an increasing function of 1/c, we face the # practical dilema of getting a solution "very close" to the ideal one (Coulomb's model) # in a reasonable amount of time. # Thus we need a small value of c nut not too small. # # The natural idea is to plot the difference between x(t) obtained with some value of c # and x(t) obtained with our "reference value" c = c__0. # # How can we do this efficiently? # (unfortunately piecewise and array/Array outputs are not allowed for parameterized solutions) # # For instance this is very slow gap := proc(c, s1, s2, n)   local s, h, k:   sol__1(parameters=[1, c]):   sol__1(s2):   h := Matrix(n, 2):   k := 1:   for s from s1 to s2 by (s2-s1)/(n-1) do     h[k, 1] := s:     h[k, 2] := eval(x(t), sol__1(s)) - eval(x(t), sol__ref(s)):     k       := k+1:   end do:   return h end proc: GAP := CodeTools:-Usage( gap(1e-5, 0, 1, 11) ): Statistics:-ScatterPlot(GAP[..,1], GAP[..,2], symbol=solidcircle, color=blue)
 memory used=399.35KiB, alloc change=0 bytes, cpu time=34.34s, real time=34.81s, gc time=0ns
 >

## Hysteresi cycle and helix movement...

Hi everyone! Currently I'm studying magnetism, and I was thinking that maybe seeing represented the helix movement of an atom with a vector v parallel to the Magnetic Field B subject to the Force of Lorentz and the Hysteresis cycle created in the magnetazation of a ferromagnetic material could help me understand it more. I tried to create the plots in Maple 2015 but I couldn't.. anyone can help me by creating those two plots?

## Problem in evaluating 2D improper integral ...

Dear Users!

I want to evaluate following integral (f[i]'s) for different values of rho as shown. But find some problem and can't find f[3], f[4], f[5] and f[6]. If the integral from direct way is hard to calculate please guide me how I can evaluate it using numerical integration?

restart; HAM := [1, 2, 3, 4, 5, 6]; U := y^6*sin(x); alpha := 1/3;

for i while i <= nops(HAM) do

rho := op(i, HAM);

f[i] := int(int(U/((x^rho-p^rho)^alpha*(y^rho-q^rho)^alpha), q = 0 .. y), p = 0 .. x)

end do

Special request to @acer @Carl Love @Kitonum @Preben Alsholm

## How to fasten this simulation?...

Hi,

I want to simulate a decision rule (let's say a production control strategy to fix the ideas) based on the outcomes of a binomial random variable.

One stage of this simulation involves operations which correspond to the NaiveProc described below. This procedure is a straightforward and rather naive translation in Maple of these operations.
It is given here to help you to understand what is to be done.

This coding being extremely slow, I implemented the required operations differently.
Given the fact that I'm only concerned by binomial samplings with outcomes 0 and 1, I basically replaced these samplings by three operations  U:=... , Got01:=... , Got0:= ... .
Each new implementation is about 200 times faster than the naive one.

My question is: Can we do still better?
For information the same simulation realized with R runs in less than 15s when REP = 10^8 (more than 100 times faster than what I'm capable to do with Maple)

 > restart:
 > with(Statistics):
 > # This is the naive coding of the problem # # It has the advantage to explain clearly (at least I hope so) what is to be done. NaiveProc := proc(REP, a, b, N)   local roll, rep, p, K, N0, N1:   roll := rand(a .. b);   N0   := 0:   N1   := 0:   for rep from 1 to REP do     p := roll();     K := Sample(Binomial(N, p), 1)[1]:     if K=0 then       N0 := N0+1:     elif K=1 then       N1 := N1+1:     end if:   end do:   N0, N1: end proc:      CodeTools:-Usage(NaiveProc(10^3, 0, 5e-4, 25000));
 memory used=78.92MiB, alloc change=68.00MiB, cpu time=1.21s, real time=6.30s, gc time=24.80ms
 (1)
 > # As the previous code is very slow I implemented several variants based on the fact # that I care only on values of K equal to 0 and 1 and that it is then useless to sample # a Binomial random variable to do this. # # Each of them is much faster than the naive coding... but can we even faster? # For comparison, the same simulation realized in R takes less than 15s to run with REP=10^8 pi := (n, p, k) -> binomial(n, k) * (1-p)^(n-k) * p^k; f1 := proc(REP, a, b, N)   local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:   Q     := Sample(Uniform(a, b), REP)^+:   P0    := pi~(N, Q, 0):   P1    := pi~(N, Q, 1):   P01   := P0 + P1:   U     := Sample(Uniform(0, 1), REP)^+:   Got01 := Select(t -> is(t <= 0), U-P01):   Got0  := Select(t -> is(t <= 0), U-P0 );   N01   := numelems(Got01):   N0    := numelems(Got0):   N1    := N01-N0:   N0, N1; end proc: CodeTools:-Usage(f1(10^5, 0, 5e-4, 25000)); print(): f2 := proc(REP, a, b, N)   local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:   Q     := Sample(Uniform(a, b), REP)^+:   P0    := pi~(N, Q, 0):   P1    := pi~(N, Q, 1):   P01   := P0 + P1:   U     := Sample(Uniform(0, 1), REP)^+:   Got01 := select[flatten](`<=`, U-P01, 0):   Got0  := select[flatten](`<=`, U-P0 , 0):   N01   := numelems(Got01):   N0    := numelems(Got0):   N1    := N01-N0:   N0, N1; end proc: CodeTools:-Usage(f2(10^5, 0, 5e-4, 25000)); print(): f3 := proc(REP, a, b, N)   local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:   Q     := Sample(Uniform(a, b), REP)^+:   P0    := pi~(N, Q, 0):   P1    := pi~(N, Q, 1):   P01   := P0 + P1:   U     := Sample(Uniform(0, 1), REP)^+:   Got01 := Vector(REP, i -> `if`(U[i] <= P01[i], 1, 0)):   Got0  := Vector(REP, i -> `if`(U[i] <= P0 [i], 1, 0)):   N01   := add(Got01):   N0    := add(Got0):   N1    := N01-N0:   N0, N1; end proc: CodeTools:-Usage(f3(10^5, 0, 5e-4, 25000));
 memory used=376.83MiB, alloc change=296.77MiB, cpu time=4.12s, real time=4.03s, gc time=430.38ms
 memory used=173.26MiB, alloc change=0 bytes, cpu time=2.30s, real time=1.98s, gc time=461.22ms
 memory used=154.19MiB, alloc change=0 bytes, cpu time=2.28s, real time=2.01s, gc time=432.40ms
 (2)
 >
 >

## Problem in execution of summation symbol ...

Dear Users!

Hope you would be fine. I have some problem in execution the last loops (highlighted as red) where sumation is present. When NN>3 it takes alot of time more than 12 hours. Is there any alternative command to reduce the query. I am waiting for your response. Thanks in advance.

restart; with(LinearAlgebra); Digits := 30; NN := 2; nu := 1; M1 := NN; M2 := NN; M3 := NN;

for k1 from 0 while k1 <= M1-1 do for k2 from 0 while k2 <= M2-1 do for k3 from 0 while k3 <= M3-1 do

SGP[M3*(M2*k1+k2)+k3+1] := simplify(sum((-1)^(k1-i1)*GAMMA(k1+i1+2*nu)*x^i1*(sum((-1)^(k2-i2)*GAMMA(k2+i2+2*nu)*y^i2*(sum((-1)^(k3-i3)*GAMMA(k3+i3+2*nu)*z^i3/(GAMMA(i3+nu+1/2)*factorial(k3-i3)*factorial(i3)), i3 = 0 .. k3))/(GAMMA(i2+nu+1/2)*factorial(k2-i2)*factorial(i2)), i2 = 0 .. k2))/(GAMMA(i1+nu+1/2)*factorial(k1-i1)*factorial(i1)), i1 = 0 .. k1)) end do end do end do;

SGPxyz := `<,>`(seq(seq(seq(SGP[M3*(M2*(i-1)+j-1)+k], k = 1 .. M3), j = 1 .. M2), i = 1 .. M1));

Lambda := `<,>`(seq(seq(seq(chi[M3*(M2*(i-1)+j-1)+k], k = 1 .. M3), j = 1 .. M2), i = 1 .. M1));

for i while i <= NN^3 do for j while j <= NN^3 do for k while k <= NN^3 do

q[i, j, k] := int(int(int(SGP[i]*SGP[j]*SGP[k]*(-x^2+x)^(nu-1/2)*(-y^2+y)^(nu-1/2)*(-z^2+z)^(nu-1/2), z = 0 .. 1), y = 0 .. 1), x = 0 .. 1) end do end do end do;

U := Matrix(NN^3, NN^3, 0);

for j while j <= NN^3 do for k while k <= NN^3 do U[j, k] := simplify(sum(chi[i1]*q[i1, j, k], i1 = 1 .. NN^3)) end do end do;

F := simplify(evalm(U));

Special request to @acer @Carl Love @Kitonum @Preben Alsholm

## A very strange result...

Hi,

I do not understand why this simple procedure evaluates so differently depending on the type of its second parameter?

A typo somewhere or a bug?

 > restart:
 > interface(version)
 (1)
 > KL := (a, b) -> (1/4)*(2*ln(a+b)*a^2+4*ln(a+b)*b*a+2*ln(a+b)*b^2-2*ln(b)*b^2-a^2-2*a*b)/a
 (2)
 > evalf(KL(1e-10, 1/2))
 (3)
 > evalf(KL(1e-10, 0.5))
 (4)
 > evalf(KL(1e-10, convert(0.5, rational)))
 (5)
 > limit(KL(a, 1/2), a=0, left); limit(KL(a, 1/2), a=0, right)
 (6)
 > limit(KL(a, 0.5), a=0, left); limit(KL(a, 0.5), a=0, right)
 (7)
 >

## Integral of a sum...

I have an integral of a sum

`Int(Sum(-(Nb*t-n*t__rev)*exp((1/2)*(L+sqrt(-4*C*L*R^2+L^2))*t/(L*R*C)-(1/2)*(Nb*t-n*t__rev)^2/(Nb^2*sigma__b^2)), n = 0 .. Nb-1), t);`

that I want to convert into a sum of integrals since Maple will not integrate the sum but it will integrate the summands. I tried IntegrationTools, but none of these tools seem to do it. It is a part of a larger expression so doing this by hand is not a real option.

Mac Dude.

## From Fortran to Maple...

Hi!

The algorithm in this PDF is written in Fortran, but unfortunately I do not known this programming languaje (actually, I am not an "expert" in progamation).

Algoritmo_Fortram.pdf  [removed by moderator. © Institute of Mathematics AS CR, 1980]

https://dml.cz/handle/10338.dmlcz/103855

Could someone please write this algorithm in Maple? Or, at least, indicate me how to do.

## Evaluation of Double Integral...

Hello Everyone
My maple is evaluating the following provided in figure. I have  attached MAPLE file. Kindly help me to solve this. Thanks in advance.

question.mw

## Is "numpoints" an effective option of SurfaceOfRe...

Hi,

While "numpoints=..." is said to be an admissible option of SurfaceOfRevolution (or VolumeOfRevolution, default value set to 50), the graphical results doesn't account for the value of numpoints.

Is this a known issue?

TIA

Surface0fRevolution.mw

 First 7 8 9 10 11 12 13 Last Page 9 of 72
﻿