vv

6875 Reputation

18 Badges

4 years, 148 days

MaplePrimes Activity


These are answers submitted by vv

restart;
PP:=proc(n::posint, kmin::posint, L::list(prime), nSolMax::posint)
# Prime partitions with elements >= kmin, starting with L
local Sol:=table(), nSol:=0, E;

E:=proc(L)
local p,sL:=add(L);
if L<>[] then p:=L[-1] else p:=nextprime(kmin-1) fi;
while p+sL<=n and nSol<=nSolMax do
  if p+sL=n and nSol<nSolMax then Sol[++nSol]:=[L[],p] fi;
  if nSol>=nSolMax then return fi;  
  E([L[],p]);
  p:=nextprime(p);
od;
end:

if L<>[] and add(L)=n then  Sol[++nSol]:=L fi;   
E(L);
entries(Sol,nolist);
end:

SingPrime:=proc(p::prime,n::posint)  
# there exist only 1 PP starting with p
evalb(nops([PP(n,p,[p],2)])=1)
end:

Q := proc(n)    # list of SingPrime's
local p, P:=select(isprime, [seq(1..n)]);
select(SingPrime, P, n)
end:

SingPrime(3,8), SingPrime(5,8);
                          true, false
for n from 2 to 500 do
    if Q(n)=[] then printf("%d ",n) fi
od;

63 161 195 235 253 425

(The procedures are not optimized.)

An implementation of Lehmer's algorithm.

Lpi := proc(x)  # Lehmer pi
  local phi, Lpi1, M:=10^7 ;
    
  phi := proc(x, a)
  local t;
  if a <= 1 then  return iquo(x + 1, 2) fi;
  t := thisproc(x, a-1) - thisproc(iquo(x, ithprime(a)), a-1);
  procname(x, a) := t;
  t;
  end;

  Lpi1 := proc(x)
  local a,b,c,s,w,i,j,L, xx:=trunc(x);
  if (trunc(x) < M) then return NumberTheory:-pi(x) fi;
  a := thisproc(iroot(xx, 4));;
  b := thisproc(isqrt(xx));
  c := thisproc(iroot(xx, 3));
  s := phi(xx,a) + (b+a-2) * (b-a+1) / 2;
  for i from a+1 to b do
    w := iquo(xx , ithprime(i));
    L := thisproc(isqrt(w));
    s := s - thisproc(w);
    if (i <= c) then
      for j from i to L do
        s := s - thisproc(iquo(w , ithprime(j))) + j - 1
      od
    fi
  od;
  s
  end;

Lpi1(x)
end:

CodeTools:-Usage(Lpi(6469693230));
memory used=6.11GiB, alloc change=80.00MiB, cpu time=55.19s, real time=50.11s, gc time=10.83s
                           300369796


 

eq := Z(x) =  3/2-9/2*exp(-2*x)-9/2*exp(-x)+1/2*int(exp(-y)*Z(x-y), y= 0..ln(2));

Z(x) = 3/2-(9/2)*exp(-2*x)-(9/2)*exp(-x)+(1/2)*(int(exp(-y)*Z(x-y), y = 0 .. ln(2)))

(1)

ZZ:=x -> A*exp(-2*x)+B*exp(-x)+C

proc (x) options operator, arrow; A*exp(-2*x)+B*exp(-x)+C end proc

(2)

expand(eval( (lhs-rhs)(eq), Z=ZZ));

(1/2)*A/(exp(x))^2+B/exp(x)+(3/4)*C-3/2+(9/2)/(exp(x))^2+(9/2)/exp(x)-(1/2)*B*ln(2)/exp(x)

(3)

coeffs(%, [exp(x),exp(2*x)]);

(3/4)*C-3/2, (1/2)*A+9/2, B+9/2-(1/2)*B*ln(2)

(4)

ABC:=solve([%],{A,B,C});

{A = -9, B = 9/(ln(2)-2), C = 2}

(5)

Zsol:=eval(ZZ(x),ABC);

-9*exp(-2*x)+9*exp(-x)/(ln(2)-2)+2

(6)

simplify(eval( (lhs-rhs)(eq), Z=unapply(Zsol,x))); # Check

0

(7)

#################################################

 

A strange fact is that intsolve finds a wrong solution but very close to it.

 

eq := Z(x) =  3/2-9/2*exp(-2*x)-9/2*exp(-x)+1/2*int(exp(-y)*Z(x-y), y= 0..ln(2));

Z(x) = 3/2-(9/2)*exp(-2*x)-(9/2)*exp(-x)+(1/2)*(int(exp(-y)*Z(x-y), y = 0 .. ln(2)))

(8)

intsolve(eq, Z(x)); # ?

Error, (in intsolve) integral equation is not linear.

 

EQ:=IntegrationTools:-Change(eq, x-y=t, t);

Z(x) = 3/2-(9/2)*exp(-2*x)-(9/2)*exp(-x)+(1/2)*(int(exp(t-x)*Z(t), t = x-ln(2) .. x))

(9)

intsolve(EQ, Z(x));

Z(x) = 3/2-(9/2)*exp(-2*x)-(9/2)*exp(-x)+(3/4)*(6*ln(2)*exp(x)-exp(2*x)+6)*exp(-x)/((ln(2)-2)*exp(x))

(10)

WrongSol:=rhs(%);

3/2-(9/2)*exp(-2*x)-(9/2)*exp(-x)+(3/4)*(6*ln(2)*exp(x)-exp(2*x)+6)*exp(-x)/((ln(2)-2)*exp(x))

(11)

plot([Zsol,WrongSol], x=0..5);

 

evalf( eval(Zsol - WrongSol,  x=1/10) );

-.938979396

(12)

evalf( eval(Zsol - WrongSol,  x=2) );

-0.932503746e-1

(13)

 


 

Download IntEq-vv.mw

Use add instead of sum i.e.

c := k -> piecewise(k = 0,b(0)^N,1 <= k,1/k/b(0)*add((q*(N+1)-k)*b(q)*c(k-q),q = 1.. k));

 

restart: 
with(GraphTheory):
L:=ImportGraph("yourpath/graph5c.g6", graph6, output=list):
DrawGraph(L);

 

The integrals are highly oscillating. t1 is convergent but not absolutely (so it does not exist as a Lebesgue integral).
Maple needs help from the user for such integrals.
I doubt Maple V was able to calculate them directly (and correctly).

Here is a special method to compute t1 (t2 goes similarly).

restart;
f1:= unapply( eval((1-exp(-(d-y)*I*x))/I/x*exp(y*I*x/(1-I*x)),  [d=3, y=9/10]), x):
F1:= u -> f1(1/u)/u^2:
h := t -> t+I*t*(1-t): # special change of variables
G1:=unapply(simplify(F1(h(t))*D(h)(t)), t):
t1:=evalf[25](2*Int(Re@f1, 0..1) + 2*Int(Re@G1, 0 .. 1));

       t1 := 3.418878800322797147227936

Edit: added a forgotten term.

Seems to be an old bug.

with(GraphTheory):
G:=Graph( Trail(1,2,3,4,5) ):
DrawGraph(G): # ok
DrawPlanar(G); # error

Probably the graph G being "collinear", it's "too planar"  :-)

 

Actually you do not have integral equations, but expressions containing integrals depending on the unknowns.
You should proceed as in the next example:

restart;
Digits:=15:
eq1 := proc(x,y) local t; int(sin(x*t)/(t^2+7), t=0..(x+y), numeric)  - 1/10. end:
eq2 := proc(x,y) local t; int(exp(x*t)/(t^2+3), t=0..(x-y), numeric) + 0.13 end:
solxy:=fsolve([eq1,eq2], [1 .. 2.5,  1.5 .. 4]);

          solxy := [1.474881121279997081543, 2.084760793500574599373]
eq1(solxy[]), eq2(solxy[]); # check
        
-.4e-15, .1e-14

restart;

alias(p = p(x, t), q = q(x, t), u = u(x, t));

p, q, u

(1)

p_x := diff(p,x) = a*p + u*q; #   (1)

diff(p, x) = a*p+u*q

(2)

q_x := diff(q,x) =  -conjugate(u)*p - a*q; # (2)     # where a is complex parameter, p_{x} mean derivative of p w.r.t x

diff(q, x) = -conjugate(u)*p-a*q

(3)

`diff/conjugate` := (z, x) ->  conjugate(diff(z,x));  # needed, x :: real

proc (z, x) options operator, arrow; conjugate(diff(z, x)) end proc

(4)

#Define, u=p^2+conjugate(q)^2,  (3)  
#Now take the derivative of (3) w.r.t x and by using (1) and (2), we get

u := p^2+conjugate(q)^2;

p^2+conjugate(q)^2

(5)

ux:=diff(u,x);

2*p*(diff(p, x))+2*conjugate(q)*conjugate(diff(q, x))

(6)

eval(ux, [p_x, q_x]);

2*p*(a*p+(p^2+conjugate(q)^2)*q)+2*conjugate(q)*conjugate(-(q^2+conjugate(p)^2)*p-a*q)

(7)

ux:=expand(%);

2*a*p^2+2*q*p^3+2*p*q*conjugate(q)^2-2*conjugate(q)^3*conjugate(p)-2*conjugate(q)*p^2*conjugate(p)-2*conjugate(q)^2*conjugate(a)

(8)

# check wrt your result
u__x := 2*(a*p^2 - conjugate(a)*conjugate(q)^2) + 2*(p^2 + conjugate(q)^2)*(p*q - conjugate(p)*conjugate(q));
simplify(ux - u__x);

2*a*p^2-2*conjugate(q)^2*conjugate(a)+2*(p^2+conjugate(q)^2)*(p*q-conjugate(p)*conjugate(q))

 

0

(9)


Download diff_conj.mw

plots:-display(seq(seq(plottools:-polygon([[m,n],[m+1/2,n+1/2],[m+1,n],[m+1/2,n-1/2]],
               color=`if`(m::odd,red,blue)), m=0..3),n=0..3), axes=none);

Note that ignoring the color, another possibility is

plots:-inequal(abs(x-floor(x)-1/2) + abs(y-floor(y)-1/2)<=1/2, x=0..4, y=0..4);

Yes, ChainTools  and also SolveTools:-SemiAlgebraic  are too slow for this problem.
But after inspection and some numeric experiments it seems that the system has solution iff:

( a=1  and  1 <= b <= sqrt(4/3) )    or    (a=b>1).

In general, such a brute-force solution is used only if a direct one is not available. This is not the case.
For t=0 the number of the roots of the equation is either <= 4 or infinity (so 6 is out of the question).

It remains t=1 which was solved in your previous post [or, just form the 4 quadratic equations (obtained from abs(u) = +- u) and find the conditions to have integer roots]. The case t=2 is similar.

I see two problems in your worksheet.

1. The distribution corresponding to the RV MIXTURE is not correctly defined. It should be:

pdf := (p,x) -> p*PDF(U1, x)*Heaviside(x+3) + (1-p)*PDF(U2, x)*Heaviside(x-1);
cdf := (p,x) -> int(pdf(p,y), y=-infinity..x);

MIXTURE := p -> RandomVariable( Distribution( PDF = (x -> pdf(p,x)), CDF = (x -> cdf(p,x)), Conditions=[p >=0, p <=1] )):

2. It seems that Sample works only for distributions for which the support is an interval and the PDF is differentiable here.
The support of  MIXTURE is (-3,-1) union (1,3).
[I don't know how to declare such a support, or even if it's possible].
So,  Sample(MIXTURE(1/3), 10);  will fail, but

Sample(MIXTURE(1/3), 10, method=[envelope, range=1..3]);  # works.

Probably it's possible to use somehow method=custom.

I have set interface(typesetting=standard), because typesetting=extended  exhibits (here too) its weakness, not related to our problem.

1. Let's take the following version of your example:

restart;
interface(typesetting=standard):
#interface(typesetting=extended):
j0:=1; n0:=0;
V:=Vector(n0):
for j from 1 to 5 do
  try
    print('j'=j);
    V(j):=j/(j-j0)
  catch:
    print("Ante",lastexception); 
    V(j):=17; # clears lastexception apparently
    print("Post",lastexception);
  end try
end do:
V, lastexception;

It is equivalent to the original one.
V(j):=17  will clear lastexception iff j0 > n0. So, it seems that the extension of a rtable is inside of some try ... end try and lastexception is reset, but only locally.
Outside try .. end try,  lastexception has always the correct value.


2. For the second example, a workaround would be:

restart;
interface(typesetting=standard);
lastexception;
5/0;
lastexception;

5/5;
lastexception;
x:=47;
lastexception;
table();
lastexception;
[1,2,3];
lastexception;
# rtable();
try  rtable(); catch : end try ;
lastexception; # Now NOT cleared!

 

The arithmetic with infinity is very delicate.
The automatic simplification must be lite, so there are compromises.
When an expression is suspected to contain infinity,  some kind of simplification is recommended.
It is strange that even is must be sometimes preceded by simplify:

eq := Pi*infinity + I*infinity = infinity + I*infinity:   #  sqrt(2) too
is(eq), is(simplify(eq));

                          false, true

(Maybe this is a consequence of the fact that a constant for complex infinity is missing in Maple.)

Edit. I add some other stange results about arithmetic with oo:

limit( (x^2+1) + x*I, x=infinity );
                            infinity
limit( (x^2+1) + x^2*I, x=infinity );
                        (1 + I) infinity
limit( exp(x) + exp(x+1)*I, x=infinity );
                   -infinity I (-exp(1) + I)
simplify(%);
                   -infinity I (-exp(1) + I)
is(%=infinity+I*infinity);
                             false

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