sand15

1576 Reputation

15 Badges

11 years, 121 days

MaplePrimes Activity


These are replies submitted by sand15

@callumneily 

See Wiki and reference herein.

 

restart

with(Statistics):

B := 2*RandomVariable(Bernoulli(1/2))-1

2*_R-1

(1)

PiEstimator := proc(Length, Rep, {path::boolean:=false})
  local samples   := Sample(B, [Length, Rep]):
  local history   := Vector[row](Length, 1) . samples;
  local estimator := 2*Length / Mean(abs~(history))^2;

  if path then
    estimator, dataplot([seq(CumulativeSum(samples[..,i]), i=1..Rep)], style=line)
  else
    estimator
  end if:
end proc:
 

Nest := 100:

PiEstimations := Vector(Nest, n -> PiEstimator(100, 1000)):

Histogram(PiEstimations);

 

# A more accute estimation

PiEstimator(10^3, 10^4)

HFloat(3.104257820893959)

(2)

# A few path


pi, pa := PiEstimator(1000, 5, path=true):
plots:-display(pa, title=typeset(Pi=evalf[4](pi)), size=[1000, 400]):

 


 

Download Pi_and_random_process.mw

 

@callumneily

 

restart

with(Statistics):

Start with an half slat of width 2b and a needle of length 2a ≤ 2b.

The kth slat covers the domain  Wk = { (x, y), x 2 (-∞, +∞) ,  y 2 [k, k+2b] }.

Consider the strip W-  = { (x, y), x 2 (-∞, +∞),  y 2 w = [-b, 0] }.
An horizontal needle is placed in W- at height Y where Y is a uniform random variable with support w.  
The needle is now rotated by a random angle T where T is a uniform random variable with support [0, p/2).  

The probability that the rotated needle intersects the line { (x, 0), x 2 (-∞, +∞) } is simply the probability
that Y + a sin(T) exceeds 0

We will first consider the case where a = b

# Without restriction:

a := 1;
b := 1;

1

 

1

(1)

Y := RandomVariable(Uniform(-b, 0)):
T := RandomVariable(Uniform(0, Pi/2)):

A := Y + a*sin(T):

Probability(A > 0)

2/Pi

(2)

f__A := PDF(A, z);

supp__A := Support(A, output=range);

N := 10^5:
E := Sample(A, N):

plots:-display(
  Histogram(E, minbins=100, style=polygon, transparency=0.7)
  ,
  plot(f__A, z=supp__A, color=blue, thickness=3)
)

`#msub(mi("f"),mi("A"))` := piecewise(z <= -1, 0, z <= 0, 2*arcsin(1+z)/Pi, z <= 1, (Pi-2*arcsin(z))/Pi, 1 < z, 0)

 

-1 .. 1

 

 


Let us consider now the case consider the case a < b = 1.

Maple (2015) can still get the solution. Indeed:

a := 'a':
assume(a > 0, a < b)
 

A := Y + a*sin(T):


Probability(A > 0);
 

2*a/Pi

(3)


But Maple (2015) fails to return the result in the general case where neither a nor b are numeric:

a := 'a':
b := 'b':
assume(b > 0):
assume(a > 0, a < b):
 

Y := RandomVariable(Uniform(-b, 0)):

A := Y + a*sin(T):


timelimit(30, Probability(A > 0));

Error, (in sprintf) time expired

 


One could get the result by explicitely computing the CDF of A derivating its expression
to get PDF(A).

Nevertheless this is a little bit lengthy and a far clever way to get the solution relies upon
the following argument. Within the strip W-  = { (x, y), x 2 (-∞, +∞),  y 2 w = [-b, 0] },
only needles whose center belong to the substrip ?D  = { (x, y), x 2 (-∞, +∞),  y 2 [-a, 0] }
can intercept the line { (x, 0), x 2 (-∞, +∞) }.
The probability Prob(Y 2 [-a, 0]) that Y takes a value in the range [-a, 0] is obviously a/b.
By probability composition rule the  probability that the needle intersects the line { (x, 0), x 2 (-∞, +∞) } 
is the product of  Prob(Y 2 [-a, 0]) by the probability that a needle with center in the substrip
 D  intersects this line.

So the result
 

Prob('intersection') = ``(2/Pi)*``('a'/'b')

Prob(intersection) = ``(2/Pi)*``(a/b)

(4)

 


 

Download BuffonNeedle.mw

(Circa 1910)

The number N of needles which intersect the edge of the slats is proportional to the needle length a and inversely proportional to the width b of the slats. So it can be written as N = P×a/b where P is some proportionality constant of which the value is to be determined.

Instead of considering straight neeedles, Borel considers circular needles of diameter d. Their perimeters, or lengths, are equal to 𝜋×d. So, with the notation above a = 𝜋×d.
Then Borel choses the particular value d = b which leads to N = P×𝜋×b/b ⟹ N = P×𝜋.

But any circle of diameter d either cuts a given slat edge twice or is tangent to two consecutive edges... which is an event of null probability.

Thus 2 = P×𝜋 ⟹ P = 2/𝜋.

@vv 

I have been taught this strict French formalism and even though I stopped doing pure math a long time ago, it has shaped my way of thinking—as it probably has for everyone who went through that phase during their math studies.

My wife, who is a little younger than I am, teaches mathematics to undergraduate and graduate students, and I asked her (before writing my previous response) whether notations such as

int(f(x), x=0..+infinity) = infinity

were acceptable or not.
She answered me that it wasn't among colleagues about her age but most common amid the "new" teacher generation, which is a constant point of discord when it is time to build the teaching programs. Kind of quarrel of the ancients and the moderns.

I was taught about coninuity and limit, thr epsilon-delts stuff when I was 15 (9th/10th grade in US educational system). Today those concepts are introduced in the second year of the mathematics and physics bachelor's program arguing that "we are going to lose to much pupils if we talk to them about continuity of a function"... and so, for the most part math student 19/20 the idea of continuity resumes to the fact the pen always slides on the sheet of paper when you draw f(x)).
I don't think the French educational math system is an exception. A few times ago there was a post about some highly ranked Maplesoft guy who had done a talk about the impoverishing of maths teaching in the US.

How can a newly appointed high school math teacher, who only discovered abstract concepts like continuity and limits (to name just a few) two years ago, possess the rigor required to teach these subjects... given that most of them are no longer even part of the course curriculum they are expected to teach?
My older daughter (31) got two master's degree in mathematics before deciding to teach maths in high school (junior and senior years). Even though she’s of course much younger than I am, she was part of the last cohorts of high school students who were taught about the epsilon-delta stuff. Sometimes she feels disheartened to see how the math curriculum has deteriorated over the past fifteen years.

Whiile some teachers lament this situation (regardless of their age, even if there are more of them among the older ones),  others simply accept it, saying, when referring to their students,“they're rubbish  and wouldn't understand,” still others (our system is almost entirely public, and teachers are civil servants paid by the state) believe it is not their role to argue against the official national curriculum (whose sole purpose is to ensure that the percentage of graduates is even higher than the previous year) finally, some have never even encountered certain basic mathematical concepts (basic at least for a 60-year-old) during their own education.

So you are right saying that "the French mathematical tradition as a whole is much more relaxed and I have seen f(oo), f(-oo) and other shorthands in many places!".
But how could it be other than that given what I wrote above.
Not making an effort to teach mathematical rigor (whether out of laziness or because that’s not what the government expects of you) inevitably leads to neglecting that very rigor and writing things like what you mentioned.

I think tools like Maple should be the last strongholds of this rigor.

@nm  @dharr  @Kitonum  

"should the underlying CAS system handle this itself and do the right thing (which would be to return undefined in this example),"

Yes, it should.
I'm French, and the French mathematics school considers writting things like 

eval(expr,x=infinity)

is a pure nonsense. Writting something like this in a maths control gives you a 0 to the question.
The correct writting in pure maths is

limit(expr, x=+infinity)


The English point of view seems far less pure and notations where the infinity symbol is consider as a number are quite common.
I have no doubt that any English mathemetician is aware that this is an abusive condensed notation. But the problem is that when a CAS accepts it, this opens the door to wrong interpretation, wrong results ... or questions on this site.
So my own opinion is that Maple should not accept people writting eval(expr,x=infinity) (or alike) and fire an error.

Another example

# Even this is incorrect (while commonly used)

int(x, x=0..+infinity)
                            infinity

# and should be written instead

limit(int(x, x=0..L), L=+infinity)
                            infinity


Note that this same difference of points of view between the English and French (maybe others) mathematics schools) appear in the notion of "domain of a function" (see Wiki for instance).
Its French definition makes a distinction between a function and an application, and is you consider an application as a function or not. Distinctions the English definition doesn't make.

Look at this simpler example: a variety of dimension 2 (the area below the curve) has a finite value while being upper bounded by a variety of lesser dimension (the length of the curve) with infinite value:

restart
f := x -> exp(-x):
int(f(x), x=0..+infinity)
                               1

L := unapply(sqrt(1+diff(f(x), x)^2), x):
int(L(x), x=0..+infinity);  # intuitively obvious result

                            infinity

I say this result is intuitively obvious because the length of the curve is necessarily larger than its projection along the x-axis... which is infinite.

This is the same as in your problem: a variety of dimension 2 (the area of gabriel's horn) of infinite 2D measure biunds a variety of dimension 3 (the volume of this horn) of finite 3D measure.
In fact the infinities you refer too are "not the same" and you cannot compare 2D and 3D measures. For instance the 3D measure of Gabriel's horn surface is 0: Would you say a finite volume is bounded with a surface of null measure?

Another point of view: the space-filling Peano curve is a curve of infinite length which is dense in the unit square. This means that a variety of dimension d may have an infinite length (wrt a d-measure) to be dense within a variety of dimension d' > d of finite d'-measure.

@acer 

"is still just as easily directly possible even if one also uses the Horner form" ... really?

What is the representation which leads to the quickest answer to the questions 
            What is the coefficient of
x6?
            What is the coefficient of x8?




Another example: What is the coefficient of x6?



But I had forgotten that you are never wrong, so I'm likely mistaken and the Horner form is surely  more suitable.

@JoyDivisionMan

Here is an alternative way to what I named the "pointwise construction" of the self-convolution of function f.
Based on the SignalProcessing:-Convolution function this computation should be much faster than the  "pointwise construction" I presented in my answer (generally numeric convolution and correlation use the FFT algorithm), but requires some carefulness to correctly rescale the result.

Here is the code (X is the same 'X' than in my reply with step=0.01).

a := Array(f~(X), 'datatype' = 'float'[ 8 ] ):
b := a:
c := SignalProcessing:-Convolution(a, b):

N := numelems(c):
K := 2*(max-min)(X)/(N-1):

#----------------------- raw convolution
# dataplot(c, style=line);

plots:-display(
  #----------------------- rescaling to get 'c' in true values
  plottools:-transform(
    (x, y) -> [x*K, y*K])
    (dataplot(c, style=line, color=red, legend="SignalProcessing:-Convolution")
  )
  ,
  plot([seq([x, Exact(x)], x in X)], legend="exact", color=blue)
  , view=[0..4, default]
)

Result

@JoyDivisionMan 

Detailed explanation using a toy problem

restart

f := x -> piecewise(x <= 0, 0, x > 2, 0, x^2);

proc (x) options operator, arrow; piecewise(x <= 0, 0, 2 < x, 0, x^2) end proc

(1)

SelfConvolution := x -> Int(f(tau)*f(x-tau), tau=0..4);

proc (x) options operator, arrow; Int(f(tau)*f(x-tau), tau = 0 .. 4) end proc

(2)


Pointwise approximation of  SelfConvolution(x)

`Approximation at point x[i]` = step * Sum('f'(h[j])*'f'(x[i]-h[j]), j=1..J);

`where:`, step =h[i+1]-h[i]

`Approximation at point x[i]` = step*(Sum(f(h[j])*f(x[i]-h[j]), j = 1 .. J))

 

`where:`, step = h[i+1]-h[i]

(3)

# A small size approximation (J=5) :

step := 1:
h    := [$(-2..2)];  

[-2, -1, 0, 1, 2]

(4)

# Select only a few number of values for x[i]

X := [$0..4];

[0, 1, 2, 3, 4]

(5)

# Apply f to each element in h

A := f~(h);

[0, 0, 0, 1, 4]

(6)

# Construct the sequence x[i]-h[1], ...x[i]-h[J]

x_shifted := x -~ h;  # shifts x by all elments of h

[x+2, x+1, x, x-1, x-2]

(7)

# Apply f to each element of x_shifted

A_shifted := f~(x_shifted );  

A_shifted := [piecewise(x <= -2, 0, 0 < x, 0, (x+2)^2), piecewise(x <= -1, 0, 0 < x-1, 0, (x+1)^2), piecewise(x <= 0, 0, 2 < x, 0, x^2), piecewise(x <= 1, 0, 0 < x-3, 0, (x-1)^2), piecewise(x <= 2, 0, 0 < x-4, 0, (x-2)^2)]

(8)

# Evaluate the Sum in relations (3) for any value of x

C := step * simplify(add(A *~ A_shifted));   # pointwise multiplication of A by A_shifted

C := piecewise(x <= 1, 0, x <= 2, (x-1)^2, x <= 3, 5*x^2-18*x+17, x <= 4, 4*(x-2)^2, 4 < x, 0)

(9)

# For comparison

simplify(value(SelfConvolution(x)));

piecewise(x < 0, 0, x <= 2, (1/30)*x^5, x < 4, 64/5-16*x+(16/3)*x^2-(1/30)*x^5, 4 <= x, 0)

(10)

# Define a fonction 'Approx' from x to C

Approx := unapply(C, x):

# Evaluate 'Approx' for all values x[i] in X.
# This gives a pointwise approximation of SelfConvolution(x):

Exact    = [ seq( [x, value(SelfConvolution(x))], x in X) ];
'Approx' = [ seq( [x, Approx(x)], x in X) ];  # Not a very good approximation but nevertheless an approximation.
                                              # To enhance it decrease the step value

Exact = [[0, 0], [1, 1/30], [2, 16/15], [3, 47/10], [4, 0]]

 

Approx = [[0, 0], [1, 0], [2, 1], [3, 8], [4, 16]]

(11)


Once you have understood the principle on this small example, do it again with a smaller value of  'step', for instance

step := 0.01;

0.1e-1

(12)

h := [seq(-2..2, step)]:
X := [seq(0..4, step)]:
A := f~(h):

x_shifted := x -~ h:
A_shifted := f~(x_shifted ):

C := step * simplify(add(A *~ A_shifted)):
S := unapply(C*step , x):

Approx := unapply(C, x):

Exact := unapply(simplify(value(SelfConvolution(x))), x);

plots:-display(
  plot([seq([x, Exact(x)], x in X)], legend="exact", color=blue)
  ,
  plot([seq([x, Approx(x)], x in X)], legend="approx", color=red)
)

proc (x) options operator, arrow; piecewise(x < 0, 0, x <= 2, (1/30)*x^5, x < 4, 64/5-16*x+(16/3)*x^2-(1/30)*x^5, 4 <= x, 0) end proc

 

 
 

 

Download Detailed_explanation_with_a_toy_example.mw


"This system of nonlinear equations was proposed by an artificial intelligence named Alice." means absolutely nothing out of context.

Alice's answers are nonsensical out of context and the last "Practical applications" item is one of the moset stupid thing I ever read.
By the way Alice forgot to mention this set of equations could also have a practical application in predicting 

She should have stayed  in Wonderland where she was more brilliant.

@Andiguys 

Question_1_sand15.mw and reply_1_sand15.mw give indeed different results because they do not solve the same equations.
To be clearer the equations you solve in your initial question and the comment to my answer are not the same.

Pay attention to what you did and read First_and_second_problem_compared.mw carefully.

@Andiguys 

Read carefully the attached file

NULL

restart

with(Optimization):

_local(Pi);

q_relation := q = (1/2)*(Ce*tau-Ci+Pr)/Cr;

q = (1/2)*(Ce*tau-Ci+Pr)/Cr

(1)

Pr_relation := Pr = -(1/2)*Ce*tau+(1/2)*q*t-(1/2)*Cd+(1/2)*Ci+1/2;

Pr = -(1/2)*Ce*tau+(1/2)*q*t-(1/2)*Cd+(1/2)*Ci+1/2

(2)

subs(Pr_relation, q_relation);

q = (1/2)*((1/2)*Ce*tau-(1/2)*Ci+(1/2)*q*t-(1/2)*Cd+1/2)/Cr

(3)

isolate((3), q)

q = (Ce*tau-Cd-Ci+1)/(4*Cr-t)

(4)

simplify(subs((4), Pr_relation), size)

Pr = ((2*Ce*tau+2*Cd-2*Ci-2)*Cr-t*(Ce*tau-Ci))/(-4*Cr+t)

(5)

``

map(simplify, solve({q_relation, Pr_relation}, {q, Pr}), size);

{Pr = ((2*Ce*tau+2*Cd-2*Ci-2)*Cr-t*(Ce*tau-Ci))/(-4*Cr+t), q = (-Ce*tau+Cd+Ci-1)/(-4*Cr+t)}

(6)


 

simplify((6)[1] - (5))

0 = 0

(7)

simplify((6)[2] - (4))

0 = 0

(8)
 

 

Download reply_1_sand15.mw

@nm 

Thanks for sharing your experience

Maybe a track: the mw file I wanted to upload contains debugger outputs: if I remove them the file can be uploaded.

For information: I get no error using Maple 2015, even for a little bit complex display.

With_Maple2015.mw

@acer 

You're right that AN/2 * BN/2 by (A*B)N/2  are not always equivalent and I had even thought to write "Assuming that AN/2 * BN/2 by (A*B)N/2  holds..."
But I was also guided by the expressions of the two terms the OP wanted to prove the equality, and without this assumption I was pretty sure the proof could not be established
So I thought that this equivalence goes without saying.

But again you're right, I sould have been detailed in my answer.

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