Maple 13 Questions and Posts

These are Posts and Questions associated with the product, Maple 13

A prime producing polynomial.


Observations on the trinomial n2 + n + 41.


by Matt C. Anderson


September 3, 2016


The story so far


We assume that n is an integer.  We focus our attention on the polynomial n^2 + n + 41.


Furthur, we analyze the behavior of the factorization of integers of the form


h(n) = n2 + n + 41                                          (expression 1)


where n is a non-negative integer.  It was shown by Legendre, in 1798 that if 0 ≤ n < 40 then h(n) is a prime number.


Certain patterns become evident when considering points (a,n) where


h(n) ≡ 0 mod a.                                             (expression 2)


The collection of all such point produces what we are calling a "graph of discrete divisors" due to certain self-similar features.  From experimental data we find that the integer points in this bifurcation graph lie on a collection of parabolic curves indexed by pairs of relatively prime integers.  The expression for the middle parabolas is –


p(r,c) = (c*x – r*y)2 – r*(c*x – r*y) – x + 41*r2.           (expression 3)


The restrictions are that 0<r<c and gcd(r,c) = 1 and all four of r,c,x, and y are integers.


Each such pair (r,c) yields (again determined experimentally and by observation of calculations) an integer polynomial a*z2 + b*z + c, and the quartic h(a*z2 + b*z + c) then factors non-trivially over the integers into two quadratic expressions.  We call this our "parabola conjecture".  Certain symmetries in the bifurcation graph are due to elementary relationships between pairs of co-prime integers.  For instance if m<n are co-prime integers, then there is an observable relationship between the parabola it determines that that formed from (n-m, n).


We conjecture that all composite values of h(n) arise by substituting integer values of z into h(a*z2 + b*z + c), where this quartic factors algebraically over Z for a*z2 + b*z + c a quadratic polynomial determined by a pair of relatively prime integers.  We name this our "no stray points conjecture" because all the points in the bifurcation graph appear to lie on a parabola.


We further conjecture that the minimum x-values for parabolas corresponding to (r, c) with gcd(r, c) = 1 are equal for fixed n.  Further, these minimum x-values line up at 163*c^2/4 where c = 2, 3, 4, ...  The numerical evidence seems to support this.  This is called our "parabolas line up" conjecture.


The notation gcd(r, c) used above is defined here.  The greatest common devisor of two integers is the smallest whole number that divides both of those integers.


Theorem 1 - Consider h(n) with n a non negative integer. 

h(n) never has a factor less than 41.


We prove Theorem 1 with a modular construction.  We make a residue table with all the prime factors less than 41.  The fundamental theorem of arithmetic states that any integer greater than one is either a prime number, or can be written as a unique product of prime numbers (ignoring the order).  So if h(n) never has a prime factor less than 41, then by extension it never has an integer factor less than 41.


For example, to determine that h(n) is never divisible by 2, note the first column of the residue table.  If n is even, then h(n) is odd.  Similarly, if n is odd then h(n) is also odd.  In either case, h(n) does not have factorization by 2.


Also, for divisibility by 3, there are 3 cases to check.  They are n = 0, 1, and 2 mod 3. h(0) mod 3 is 2.  h(1) mod 3 is 1. and h(2) mod 3 is 2.  Due to these three cases, h(n) is never divisible by 3.  This is the second column of the residue table.


The number 0 is first found in the residue table for the cases h(0) mod 41 and h(40) mod 41.  This means that if n is congruent to 0 mod 41 then h(n) will be divisible by 41.  Similarly, if n is congruent to 40 mod 41 then h(n) is also divisible by 41.

After the residue table, we observe a bifurcation graph which has points when h(y) mod x is divisible by x.  The points (x,y) can be seen on the bifurcation graph.


< insert residue table here >


Thus we have shown that h(n) never has a factor less than 41.


Theorem 2


Since h(a) = a^2 + a + 41, we want to show that h(a) = h( -a -1).


Proof of Theorem 2

Because h(a) = a*(a+1) + 41,

Now h(-a -1) = (-a -1)*(-a -1 +1) + 41.

So h(-a -1) = (-a -1)*(-a) +41,

And h(-a -1) = h(a).

Which was what we wanted.

End of proof of theorem 2.


Corrolary 1

Further, if h(b) mod c ≡ = then h(c –b -1) mod c ≡ 0.


We can observe interesting patterns in the “graph of discrete divisors” on a following page.



Dear Colleagues,


I am not sure if there exist a simple way to handle the issues I am facing. I am trying to obtain numeric roots for a polynomial f(x,a). I know for sure that there can be many roots depending on the value of parameter a. However, I cannot say for sure how many roots are possible for each value of parameter a. Some of these roots are complex numbers. Also, I need to choose only those roots that have following properties:

1. They are real.

2. f(x*,a) i.e., function value at a root is positive. 


How do I solve f(x,a)=0 to store all roots in a set? Furthermore, how do I select and print roots that have the properties mentioned above? Is there a way to do filtering of a set specifying properties of the members of the set? Please suggest. Your help is highly appreciated.







i want a scheme of fractional differential equation so that i solve my questions and make a code of it.

please provide me the scheme

i have an example, u[t] = u[xx]^2+u[yy]^2+u[zz]^2 with subject to b.c. u[0](x,y,z,t):=2*sin(x)*sin(y)*sin(z)
i used adomian method to solve this P.D.E, but i failed to construct a code of 2D P.D.E.
kindly help me in this regard

Typically sets are created like:



and then you can carryout A union B or B\A


what if  you wanted to create the set as

A:={values in some three dimensional space};

B:={volume, based on values taken from A};

Can these relationships be set up in Maple? If so, how? If there are commands that specifically handle these types of sets, what does maple call them?  I've seen the term 'set function' but what might Maple call them?

Note: I am not even sure i 'tagged' this correctly because it I am not sure the proper terms for these functions/sets.

Thanks in advance for any help.

Hi MaplePrimes,


I try again to ask about a procedure I wrote.  Maybee someone can spot my error.

Here is the broken procedure - 

cs := proc (n::integer)::integer;
local count;
description "Find the number of steps for a Collatz sequence to reach 1";
count := 1;
while 1 < n do
if `mod`(n, 2) = 0 then n := (1/2)*n else n := 3*n+1 end if;
count := count+1
end do;
return count
end proc




Hi everyone!

I have a problem solving the nonlinear ode (as attached below). I got this error ---> Error, (in fproc) unable to store '-1.32352941215398+(-0.441176470717993e-1, -0.)' when datatype=float[8]

1) Could someone please explain to me what does the unable store .... error means? 

and i will be grateful if you could help me finding the solution out. Thanks in advance

How can I plot the complex function f(z)=1/(1-z) for |z|<1 into maple code? (z is a complex number)


Best regards,


How can I modify the appearance of the arrowheads on the vectors displayed in phaseportrai? In particular, how can I "fill in" the arrowheads so that the arrowhead is not just an outline?

My code is:


Thank you

Hi everyone,

I'm kinda new here, and I really hope you guys can help me through this. In my new case study, after some revision, i thought i might be trying to implement a shooting method. I tried my best to make it work/understand but i couldn't get to any result.

So, as attached (i re-do PV Satya Naraya's paper first to be more understand but .....)


Here is my questions and the worksheet:

1) really stuck in mind - what is the purpose of shooting method for some related study?

2) what is the meaning of error .............'use midpoint method intead" 

3) Worksheet -

Thanks in advanced. Really hope that someone can help/teach me how to solve the boundary value problem by shooting method. 



restart; with(plots); lambda := 1.0; m := 2.0; M := 2; R := .1; Pr := .75; G := .1; Sc := .6; Kr := .2; blt := 5

Eq1 := diff(f(eta), eta, eta, eta)+(1+lambda)*(f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2)-(1+lambda)*M*(diff(f(eta), eta))+beta*((diff(f(eta), eta, eta))^2-f(eta)*(diff(f(eta), eta, eta, eta, eta))) = 0;

diff(diff(diff(f(eta), eta), eta), eta)+2.0*f(eta)*(diff(diff(f(eta), eta), eta))-2.0*(diff(f(eta), eta))^2-4.0*(diff(f(eta), eta))+beta*((diff(diff(f(eta), eta), eta))^2-f(eta)*(diff(diff(diff(diff(f(eta), eta), eta), eta), eta))) = 0



Eq2 := (1+(4/3)*R)*(diff(theta(eta), eta, eta))+Pr*(f(eta)*(diff(theta(eta), eta))-m*(diff(f(eta), eta))*theta(eta)+G*theta(eta)) = 0;

1.133333333*(diff(diff(theta(eta), eta), eta))+.75*f(eta)*(diff(theta(eta), eta))-1.500*(diff(f(eta), eta))*theta(eta)+0.75e-1*theta(eta) = 0


Eq3 := diff(phi(eta), eta, eta)+Sc*(f(eta)*(diff(phi(eta), eta))-m*(diff(f(eta), eta))*phi(eta)-Kr*phi(eta)) = 0;

diff(diff(phi(eta), eta), eta)+.6*f(eta)*(diff(phi(eta), eta))-1.20*(diff(f(eta), eta))*phi(eta)-.12*phi(eta) = 0


bcs1 := f(0) = 0, (D(f))(0) = 1, (D(f))(blt) = 0, (D(D(f)))(blt) = 0, theta(0) = 1, theta(blt) = 0, phi(0) = 1, phi(blt) = 0;

f(0) = 0, (D(f))(0) = 1, (D(f))(5) = 0, ((D@@2)(f))(5) = 0, theta(0) = 1, theta(5) = 0, phi(0) = 1, phi(5) = 0


L := [1.0, 1.5, 2.0, 2.5];

[1.0, 1.5, 2.0, 2.5]


for k to 4 do R := dsolve(eval({Eq1, Eq2, Eq3, bcs1}, beta = L[k]), [f(eta), theta(eta), phi(eta)], numeric, output = listprocedure); Y || k := rhs(R[3]); YA || k := rhs(R[6]); YB || k := rhs(R[5]); YC || k := -rhs(R[8]) end do

Error, (in dsolve/numeric/bvp) system is singular at left endpoint, use midpoint method instead









  Hi, there

How can I draw the excircles, incircles,circumcircle and their centers of a triangle simultaneously with maple13 in a geometric plot? please specify the commands.

many thanks for your help


I am using maple 13 to found Eingenvalues of an hermitian matrix :

> [lambda3+lambda4,0,0,0,0,0,lambda4/sqrt(2),0,0,I*lambda4/sqrt(2)],
> [0,lambda3/4,0,0,0,0,0,0,0,0],
> [0,0,lambda3/4,0,0,0,0,0,0,0],
> [0,0,0,lambda3/4,0,0,0,0,0,0],
> [0,0,0,0,lambda3/4,0,0,0,0,0],
> [0,0,0,0,0,lambda3,0,0,0,0],
> [lambda4/sqrt(2),0,0,0,0,0,lambda3/2,0,0,0],
> [0,0,0,0,0,0,0,lambda2,0,0],
> [0,0,0,0,0,0,0,0,lambda2,0],
> [-I*lambda4/sqrt(2),0,0,0,0,0,0,0,0,lambda4/2]
> ]);


my surprise is that maple gives me 8 correct solutions an 2 complex eigenvalues which are not acceptable (we now that the eigenvalues for an hermitian matrix are all real) .

To understand the output of maple, first,  I suspect that the complex part of the roots was null but without success I haven't found how to do it zero...

is it a bug? Thanks a lot to cooperation

Dear Maple researchers


I have a problem in solving a system of odes that resulted from discretizing, in space variable, method of lines (MOL).

The basic idea of this code is constructed from the following paper:

If kindly is possible, please tell me whas the solution of this problem.

With kin dregards,

Emran Tohidi.

My codes is here:

> restart;
> with(orthopoly);
print(`output redirected...`); # input placeholder
> N := 4; Digits := 20;
print(`output redirected...`); # input placeholder

> A := -1; B := 1; rho := 3/4;
> g1 := proc (t) options operator, arrow; 1/2+(1/2)*tanh((1/2)*(A-(2*rho-1)*t/sqrt(2))/sqrt(2)) end proc; g2 := proc (t) options operator, arrow; 1/2+(1/2)*tanh((1/2)*(B-(2*rho-1)*t/sqrt(2))/sqrt(2)) end proc;
print(`output redirected...`); # input placeholder
> f := proc (x) options operator, arrow; 1/2+(1/2)*tanh((1/2)*x/sqrt(2)) end proc;
print(`output redirected...`); # input placeholder
> uexact := proc (x, t) options operator, arrow; 1/2+(1/2)*tanh((1/2)*(x-(2*rho-1)*t/sqrt(2))/sqrt(2)) end proc;
print(`output redirected...`); # input placeholder
> basiceq := simplify(diff(uexact(x, t), `$`(t, 1))-(diff(uexact(x, t), `$`(x, 2)))+uexact(x, t)*(1-uexact(x, t))*(rho-uexact(x, t)));
print(`output redirected...`); # input placeholder
> alpha := 0; beta := 0; pol := P(N-1, alpha+1, beta+1, x); pol := unapply(pol, x); dpol := simplify(diff(pol(x), x)); dpol := unapply(dpol, x);
print(`output redirected...`); # input placeholder
> nodes := fsolve(P(N-1, alpha+1, beta+1, x));
> xx[0] := -1;
> for i to N-1 do xx[i] := nodes[i] end do;
print(`output redirected...`); # input placeholder
> xx[N] := 1;
> for k from 0 to N do h[k] := 2^(alpha+beta+1)*GAMMA(k+alpha+1)*GAMMA(k+beta+1)/((2*k+alpha+beta+1)*GAMMA(k+1)*GAMMA(k+alpha+beta+1)) end do;
print(`output redirected...`); # input placeholder
> w[0] := 2^(alpha+beta+1)*(beta+1)*GAMMA(beta+1)^2*GAMMA(N)*GAMMA(N+alpha+1)/(GAMMA(N+beta+1)*GAMMA(N+alpha+beta+2));
print(`output redirected...`); # input placeholder
> for jj to N-1 do w[jj] := 2^(alpha+beta+3)*GAMMA(N+alpha+1)*GAMMA(N+beta+1)/((1-xx[jj]^2)^2*dpol(xx[jj])^2*factorial(N-1)*GAMMA(N+alpha+beta+2)) end do;
print(`output redirected...`); # input placeholder
> w[N] := 2^(alpha+beta+1)*(alpha+1)*GAMMA(alpha+1)^2*GAMMA(N)*GAMMA(N+beta+1)/(GAMMA(N+alpha+1)*GAMMA(N+alpha+beta+2));
print(`output redirected...`); # input placeholder
> for j from 0 to N do dpoly1[j] := simplify(diff(P(j, alpha, beta, x), `$`(x, 1))); dpoly1[j] := unapply(dpoly1[j], x); dpoly2[j] := simplify(diff(P(j, alpha, beta, x), `$`(x, 2))); dpoly2[j] := unapply(dpoly2[j], x) end do;
print(`output redirected...`); # input placeholder
print(??); # input placeholder
> for n to N-1 do for i from 0 to N do BB[n, i] := sum(P(jjj, alpha, beta, xx[jjj])*dpoly2[jjj](xx[n])*w[i]/h[jjj], jjj = 0 .. N) end do end do;
> for n to N-1 do d[n] := BB[n, 0]*g1(t)+BB[n, N]*g2(t); d[n] := unapply(d[n], t) end do;
print(`output redirected...`); # input placeholder
> for nn to N-1 do F[nn] := simplify(sum(BB[nn, ii]*u[ii](t), ii = 1 .. N-1)+u[nn](t)*(1-u[nn](t))*(rho-u[nn](t))+d[nn](t)); F[nn] := unapply(F[nn], t) end do;
print(`output redirected...`); # input placeholder
> sys1 := [seq(d*u[q](t)/dt = F[q](t), q = 1 .. N-1)];
print(`output redirected...`); # input placeholder
[d u[1](t)                                                                
[--------- = 40.708333333333333334 u[1](t) + 52.190476190476190476 u[2](t)
[   dt                                                                    

                                                                  2          3
   + 39.958333333333333334 u[3](t) - 1.7500000000000000000 u[1](t)  + u[1](t)

   + 7.3392857142857142858

   - 3.6696428571428571429 tanh(0.35355339059327376220

   + 0.12500000000000000000 t) - 3.6696428571428571429 tanh(
                                                     d u[2](t)   
-0.35355339059327376220 + 0.12500000000000000000 t), --------- =
-20.416666666666666667 u[1](t) - 25.916666666666666667 u[2](t)

                                                                  2          3
   - 20.416666666666666667 u[3](t) - 1.7500000000000000000 u[2](t)  + u[2](t)

   - 3.7500000000000000000

   + 1.8750000000000000000 tanh(0.35355339059327376220

   + 0.12500000000000000000 t) + 1.8750000000000000000 tanh(
                                                     d u[3](t)                
-0.35355339059327376220 + 0.12500000000000000000 t), --------- = 29.458333333\

  333333333 u[1](t) + 38.476190476190476190 u[2](t)

                                                                  2          3
   + 30.208333333333333333 u[3](t) - 1.7500000000000000000 u[3](t)  + u[3](t)

   + 5.4107142857142857144

   - 2.7053571428571428572 tanh(0.35355339059327376220

   + 0.12500000000000000000 t) - 2.7053571428571428572 tanh(
-0.35355339059327376220 + 0.12500000000000000000 t)]
> ics := seq(u[qq](0) = evalf(f(xx[qq])), qq = 1 .. N-1);
print(`output redirected...`); # input placeholder
    u[1](0) = 0.38629570659055483825, u[2](0) = 0.50000000000000000000,

      u[3](0) = 0.61370429340944516175
> dsolve([sys1, ics], numeic);
Error, (in dsolve) invalid input: `PDEtools/sdsolve` expects its 1st argument, SYS, to be of type {set({`<>`, `=`, algebraic}), list({`<>`, `=`, algebraic})}, but received [[d*u[1](t)/dt = (20354166666666666667/500000000000000000)*u[1](t)+(13047619047619047619/250000000000000000)*u[2](t)+(19979166666666666667/500000000000000000)*u[3](t)-(7/4)*u[1](t)^2+u[1](t)^3+36696428571428571429/5000000000000000000-(36696428571428571429/10000000000000000000)*tanh(1767766952966368811/5000000000000000000+(1/8)*t)-(36696428571428571429/10000000000000000000)*tanh(-1767766952966368811/5000000000000000000+(1/8)*t), d*u[2](t)/dt = -(20416666666666666667/1000000...

How difficult is it to simulate gravitational influences and perturbing effects on celestial orbits with Maple? Could this syntax be altered without excessive changes to consider these aspects?

Are there somewhere worksheets to take a look at as an introduction and to see how such goals would be approached and implemented?

Some time ago, I have used the Virtual 3D Solar System code to plot a few interesting figures with asteroids. It can be found here: . I plotted the asteroid belt as a closed, warped surface according to how the asteroids are inclined, with the arguments of perihelion and longitude of ascending node where they cross the equatorial plane to the other celestial hemisphere. The plot was only done with axes in units of distances, for convenience's sake in astronomical units of 1 AU = 150 million km.

But I am strongly wondering if it were somehow possible to alter or adjust the code such that there is a zylinder placed around such plot where one can see celestial coordinates (Rectascension, Declination) on the inside, so that if one wants to know what coordinates an orbit of a specific asteroid has that one can look in the chart and go out on the balcony and set the telescope to these coordinates. Of course the asteroid most likely is not there on his entire orbital path, but you have at least the location where the asteroid CAN be. Would this be possible to plot?

This probably doesn't work if a closed cylinder is placed around the Virtual Solar System. But how about scrolling the cylinder and the orbits flat on a 2D plot?

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