mmcdara

3000 Reputation

16 Badges

5 years, 187 days

MaplePrimes Activity


These are answers submitted by mmcdara

The method is described in the help pages, please look to the attached file (the green texts are clickable) for an exampleDownload hyperlink.mw

Even if I'm not familiar with this package, the Physics[Setup] help page says that vectordisplay should do the job:

vectordisplay   none, arrow, bold, or halfarrow.
When mathematicalnotation is set to true, the vectordisplay setting causes the Maple Standard GUI to display non-projected vectors with an arrow or halfarrow on top, or in bold, as specified.

An example (Maple 2015)

restart:
CountPrimes := proc(p::integer, q::integer)
  local k := 0:
  local i := nextprime(p):
  while i < q do
    k := k+1:
    i := nextprime(i):
  end do:
  return k
end proc:
      
start := 1:
for p in [seq(1+1000*k, k=0..9)] do
  printf("range %4d..%5d : nb of primes = %d\n", p, p+999, CountPrimes(p, p+999))
end do:

range    1.. 1000 : nb of primes = 168
range 1001.. 2000 : nb of primes = 135
range 2001.. 3000 : nb of primes = 127
range 3001.. 4000 : nb of primes = 119
range 4001.. 5000 : nb of primes = 118
range 5001.. 6000 : nb of primes = 114
range 6001.. 7000 : nb of primes = 117
range 7001.. 8000 : nb of primes = 106
range 8001.. 9000 : nb of primes = 110
range 9001..10000 : nb of primes = 111

Concerning the plot, here is a suggestion 

res := NULL:
start := 1:
for p in [seq(1+1000*k, k=0..9)] do
  res := res, plottools:-rectangle([p, 0], [p+999, CountPrimes(p, p+999)], transparency=0.7)
end do:
plots:-display(res)

To conclude:
"Can i make a procedure what calculates the latest mersenne prime number ?"  what does "latest" mean?
Advice: look at numtheory[mersenne] (Maple 2015, think it is numbertheory[mersenne] in more recent versions)

I believe you misinterpreted the phrase "This routine returns the "normalized" form of a linear differential equation, a differential operator list,...".

As written below it a differential operator list is not a list of odes but the list of the coefficients of the successive derivatives of the dependent variable of an ode.

As written in the DEnormal help page this list can be provided directly from the user, or derived from the ode by using convertAlg.
See the code bellow which is just a a commented version of the one given in this help page

restart:
with(DEtools):

# case of an ode
DE := 21*(-2*x^3)/(3*x^2-5*x+2)*y(x) + 42*(x^2)/(x+1)*x*(x-1)*D(y)(x) +
50*x^3/(x-1)^3*(D@@2)(y)(x) = x*sin(x);
DEnormal(DE,x,y(x));


# case of an operator list (here the operator list which corresponds to DE)
DE2 := convertAlg(DE,y(x));
DEnormal(DE2,x);

In the  convertAlg help page it's said that it "returns the coefficient list form for a linear ODE" (not of system of linear odes).
You can apply convertAlg and DEnormal for each ode in L this way

L := [diff(z(t), t, t)*m = 0, diff(x(t), t, t)*m - B*diff(y(t), t)*q = 0, diff(y(t), t, t)*m + B*diff(x(t), t)*q = 0] ;
V := [z(t), x(t), y(t)];
A := [seq(convertAlg(L[i], V[i]), i=1..3)];
N := DEnormal~(A, t);

and finally obtain the normal form of each ode this way

map(i -> add(N[i][1][j]*~(D@@(j-1))(op(0, V[i]))(t), j=1..nops(N[i][1])) = N[i][2], [$1..3]):
convert~(%, diff)
     [                                       / d      \    
     [                                     B |--- y(t)| q  
     [ d  / d      \       d  / d      \     \ dt     /    
     [--- |--- z(t)| = 0, --- |--- x(t)| = --------------, 
     [ dt \ dt     /       dt \ dt     /         m         

                            / d      \  ]
                          B |--- x(t)| q]
        d  / d      \       \ dt     /  ]
       --- |--- y(t)| = - --------------]
        dt \ dt     /           m       ]


But it's much simpler to do this

`Normal form of L`:= DEnormal~(L, t, V);

DEnormal.mw

Here is an example 

restart:

f := x -> x^2;
g := x -> x^3-1;
h := x -> abs(f(x)-g(x));

u := x -> ['x', 'f(x)', 'g(x)', 'h(x)'];

X := [seq(k, k = 1 .. 5)]:
M := convert([u(x), eval(u~(X))[]], Matrix);

# or maybe
M := convert([eval(u(x)), eval(u~(X))[]], Matrix);

# M can be saved using 
#    ExcelTools:-Export
#    ExportMatrix
#    ...
#
# For instance

n := cat(currentdir(), /"MyMatrix":
ExportMatrix(n, M);



Create_rable_mmcdara.mw

You could also use series(rhs(sol), x=0, 8) instead of asympt in your second example.

For your first example here is the Maple's solution based on a series expansion of the formal solution.
You will see it differs from MMA's from many points and I wasn't capable to transform th Maple's solution into MMA's one.

series.mw

But I guess you already did this yourself?

For the inverse Laplace transform of a function f(s) might exist, f(s) must vanish to 0 faster than
|s^n * f(s)|  (n strictly positive integer) as s tend to infinity.
This is not the case of d(s)


I don't say it's more efficient than what @tomleslie did, but I do believe it is more elegant and closer yo the formal problem:

 

restart:

with(Statistics):

# Be careful:
# In the conventional parameterization of an exponential distribution
# its parameter (p) has the same dimension as the realizations of the random variable.
# Look at this

E := RandomVariable(Exponential(p)):
PDF(E, t)
 

piecewise(t < 0, 0, exp(-t/p)/p)

(1)

# with this is mind your RV should be defined this way

E__X := RandomVariable(Exponential(1/lambda)):
E__Y := RandomVariable(Exponential(1/mu)):
 

# Instead of explicitely write something like
#      assuming lambda > 0, mu > 0
# it's better (IMO) to use the conditions imposed to lambda and mu
# as parameters of exponential distributions:

att := attributes(E__X)[3]: C__X := att:-Conditions[];
att := attributes(E__Y)[3]: C__Y := att:-Conditions[];

0 < 1/lambda

 

0 < 1/mu

(2)

# Then the expected result simply writes

 Probability(E__Y > E__X)  assuming C__X, C__Y

lambda/(lambda+mu)

(3)

 


 

Download The_probabilistic_way.mw

Obviously 7, 8, 12.. are othees among an infinity
 

restart:
sols := isolve(x^2+1=5*k):
S := op~(2, [sols]);

                 [x = 3 - 5 _Z1, x = 2 - 5 _Z1]

seq(op(eval(S, _Z1=n)), n=-10..0):
sort(rhs~([%]))
[2, 3, 7, 8, 12, 13, 17, 18, 22, 23, 27, 28, 32, 33, 37, 38, 42, 

  43, 47, 48, 52, 53]

 


Look to the attached file
proposal.mw

Like @Mac Dude I'm not sure I truly understand what you want.
Maybe something like this?
Maybe_This.mw

Here is a very simple version (continuity corrections do exist but are not implemented here).
Let E some event whose outcomes are 0 and 1

  • n = sample size
  • k = numer of times E=1 is observed amid these n observations
  • p = the hypothetized probability that E=1
  • alpha = the risk of first kind
  • side = either "L" (left tail), "R", right tail or "LR" two-tailed (symmetric)
     

restart

ExactBinomialTest := proc(n, k, p, alpha, side)
  local B, beta, ev, lb, ub, q:
  uses Statistics:
  B    := RandomVariable(Binomial(n, p)):
  ev   := n*p:
  beta := `if`(side="LR", alpha/2, alpha):
  lb   := Quantile(B, beta):
  ub   := Quantile(B, 1-beta):
  #printf("The expected value of k given P=%a and N=%d is %d\n", p, n, ev):


  if k <= ev then
    q := Probability(B <= k, numeric)
  else
    q := Probability(B >= k, numeric)
  end if:
  if k < lb then
    printf("The null hypothesis P0=%a is rejected\n", p):
  elif k > ub then
    printf("The null hypothesis P0=%a is rejected\n", p):
  else
    printf("The null hypothesis P0=%a is accepted\n", p):
  end if:
  printf("p_value = %1.3f\n", q):
end proc:

ExactBinomialTest(100, 41, 1/2, 0.05, "L");
ExactBinomialTest(100, 41, 1/2, 0.05, "LR"):
ExactBinomialTest(100, 41, 1/2, 0.05, "R")

The null hypothesis P0=1/2 is rejected

p_value = 0.044
The null hypothesis P0=1/2 is accepted
p_value = 0.044
The null hypothesis P0=1/2 is rejected
p_value = 0.044

 

ExactBinomialTest(100, 40, 1/3, 0.10, "L");
ExactBinomialTest(100, 40, 1/3, 0.10, "LR"):
ExactBinomialTest(100, 40, 1/3, 0.10, "R")

The null hypothesis P0=1/3 is rejected

p_value = 0.097
The null hypothesis P0=1/3 is accepted
p_value = 0.097
The null hypothesis P0=1/3 is rejected
p_value = 0.097

 

ExactBinomialTest(10^6, 80, 1e-4, 0.10, "L")

The null hypothesis P0=.1e-3 is rejected

p_value = 0.023

 

 

 

Download ExactBinomialTest.mw

I was writting the attached file as @dharr was posting his own answer.

This is basically the same thing, probably more verbose in my case
 

restart; with(plottools); with(plots)

 

In John Stillwell's book "The Four Pillars of Geometry" he declares the following to be the equations of non-Euclidean "lines" where A,B and C are all non-zero


What Maple code will plot this equation for given values of A, B and C?

Equation := A*z*conjugate(z)+B*(z+conjugate(z))+C = 0:

# assumptions : A::real, B::real, C::real

expand(eval(Equation, z=x+I*y)) assuming x::real, y::real;

A*x^2+A*y^2+2*B*x+C = 0

(1)

expand(%/A);
Student:-Precalculus:-CompleteSquare(%, x);
eq := select(has, lhs(%), {x, y}) = -select(`not`@`has`, lhs(%), {x, y})

x^2+y^2+2*B*x/A+C/A = 0

 

(x+B/A)^2+y^2+C/A-B^2/A^2 = 0

 

(x+B/A)^2+y^2 = -C/A+B^2/A^2

(2)

local gamma:
eq := algsubs(C/A=gamma, algsubs(B/A=beta, eq));

Warning, A new binding for the name `gamma` has been created. The global instance of this name is still accessible using the :- prefix, :-`gamma`.  See ?protect for details.

 

(x+beta)^2+y^2 = beta^2-gamma

(3)

# Case rhs(eq) > 0
#
# eq is the equation of a circle of center (-beta, 0) iif gamma verifies

solve(rhs(eq) >= 0, gamma)

{gamma <= beta^2}

(4)

# Case rhs(eq) = 0
#
# the only couple (x,y) which verifies eq when rhs(eq)=0 is

x=-beta, y=0;  # as x, y and beta were assumed real

x = -beta, y = 0

(5)

# Case rhs(eq) < 0
#
# if rhs(eq) < 0, then eq is of the form (x+beta)^2+y^2 = -r^2 with r > 0

sol := solve(eq, y) assuming rhs(eq) < 0, x::real, y::real;

(-2*beta*x-x^2-gamma)^(1/2), -(-2*beta*x-x^2-gamma)^(1/2)

(6)

# -2*beta*x - x^2 - gamma = -2*beta*x - x^2 - gamma + beta^2 - beta^2
# -2*beta*x - x^2 - gamma = -(2*beta*x + x^2 + beta^2) + (beta^2 - gamma)
# -2*beta*x - x^2 - gamma = -(x + beta)^2 + (beta^2 - gamma)
# -2*beta*x - x^2 - gamma < 0  for each real x
#
# thus each member of sol is a complex number which contradicts the assumption y::real.
#
# Case rhs(eq) < 0 has no (x, y) real solution

pl := proc(a, b, c, h)
  local __beta  := b/a:
  local __gamma := c/a:
  local r := __beta^2 - __gamma;
  if r > 0 then
    implicitplot(eval(eq, [beta=__beta, gamma=__gamma]), x=-h..h, y=-h..h, grid=[4*h^2, 4*h^2], scaling=constrained)
  elif r = 0 then
    plot([[-__beta, 0]], style=point, symbol=circle, symbolsize=15, scaling=constrained)
  else
    WARNING("no real solution");
  end if:
end proc:


pl(1, 2, 1, 5);

 

pl(1, 2, 4, 5);

 

pl(1, 2, 5, 5);

Warning, no real solution

 

 


 

Download NonEuclideanLines_mmcdara.mw

 

Curiously the definitions of a vector basis seem to differ in English and in French.

In English a set F of vectors of some space E is a basis if every element of E can be written as a unique linear combination of elements of F.

In French a set ("family") F which verifies this condition is termed a "generating family".
A family F={u, v} is termed "free" if then a*u+b*v=0 ==> a=b=0.
And F is a basis iif it is a generating family and a free family. 

If the number of elements of F equals the dimension of V, then it's enough to veriffy that the family is a generating one or that it is a free one.
Generally it's easier to verify this second condition (in such a case the family is termed "maximally free"):

restart;
u := <1, 1>:
v := <1, 2>:

# Is {u,v} a free family?
cond := lambda*~u + mu*~v =~ <0, 0>:
solve(convert(cond, list), [lambda, mu])
                     [[lambda = 0, mu = 0]]

Correction is in 1D input mode
 

restart

with(LinearAlgebra); A := Matrix(5, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = tau*`&Pi;u`/mu, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = `&varpi;`*`&Pi;g`/mu}); B := Matrix(5, 5, {(1, 1) = mu, (1, 2) = 0, (1, 3) = 0, (1, 4) = gamma, (1, 5) = 0, (2, 1) = 0, (2, 2) = mu+omega, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = mu+sigma1, (3, 4) = theta, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = alpha2+gamma+mu, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = sigma1+alpha2, (5, 4) = 0, (5, 5) = theta}); C := A.(1/B); Rank(C); evs := Eigenvalues(C); eig := op(`minus`({entries(evs, nolist)}, {0}))

A := Matrix(5, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = tau*`&Pi;u`/mu, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = varpi*`&Pi;g`/mu})

 

B := Matrix(5, 5, {(1, 1) = mu, (1, 2) = 0, (1, 3) = 0, (1, 4) = gamma, (1, 5) = 0, (2, 1) = 0, (2, 2) = mu+omega, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = mu+sigma1, (3, 4) = theta, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = alpha2+gamma+mu, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = sigma1+alpha2, (5, 4) = 0, (5, 5) = theta})

 

C := Matrix(5, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = tau*`&Pi;u`/(mu*(mu+sigma1)), (3, 4) = -tau*`&Pi;u`*theta/(mu*(alpha2+gamma+mu)*(mu+sigma1)), (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = -varpi*`&Pi;g`*(sigma1+alpha2)/(mu*(mu+sigma1)*theta), (5, 4) = varpi*`&Pi;g`*(sigma1+alpha2)/(mu*(mu+sigma1)*(alpha2+gamma+mu)), (5, 5) = varpi*`&Pi;g`/(mu*theta)})

 

2

 

evs := Vector(5, {(1) = 0, (2) = 0, (3) = 0, (4) = tau*`&Pi;u`/(mu*(mu+sigma1)), (5) = varpi*`&Pi;g`/(mu*theta)})

 

tau*`&Pi;u`/(mu*(mu+sigma1)), varpi*`&Pi;g`/(mu*theta)

(1)

params := {`&Pi;g` = 2.2, `&Pi;u` = 3.4, alpha2 = .33, mu = .2041, omega = .5, sigma1 = .72, tau = .33, theta = .9, varpi = 0.96e-1};

{`&Pi;g` = 2.2, `&Pi;u` = 3.4, alpha2 = .33, mu = .2041, omega = .5, sigma1 = .72, tau = .33, theta = .9, varpi = 0.96e-1}

(2)

plot(eval(eig, remove(has, params, `&varpi;`)), `&varpi;` = 0 .. .5, filled = true, thickness = 4, color = blue, transparency = .4)

Error, (in plot) unexpected options: [11.97669988*varpi, varpi = 0 .. .5]

 

# what are you trying to plot?

eval(eig, remove(has, params, varpi))

5.948820737, 11.97669988*varpi

(3)

# two ways :
#  1/ plot these 2 quantities

plot([eval(eig, remove(has, params, varpi))], varpi = 0 .. .5, filled = true, thickness = 4, color = [blue, red], transparency = .4)

 


#  2/ plot only the last one

plot(eval(eig, remove(has, params, varpi))[2], varpi = 0 .. .5, filled = true, thickness = 4, color = blue, transparency = .4)

 

 


 

Download graph_mmcdara.mw

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