sand15

1449 Reputation

15 Badges

11 years, 80 days

MaplePrimes Activity


These are replies submitted by sand15

Error messages are pretty clear: you have to specify an approximate solution.

For more detais execute help(dsolve,numeric_bvp,advanced) and help(dsolve,numeric,BVP) in a worksheet and search for the word approxsoln.

Defining an approximate solution can be quite tricky and the best way to do this is to construct a simplified version your problem using a simplified physics (do a dimention analysis, identify second order effects, discard them), to solve it, and to use this solution as approxsoln (which can be either a procedure or an array).
Your are likely the better placed to do this because it is you who knows the physics behind your equations.

For instance, what is order of magnitude of
compared to other terms? Can you simplify it?

Or can you approximate BCs like?
Can you simplify them from a physical point of view?
As a rule a BC on the first derivative which nonlinearly depends on the second derivative makes a BC problem complicated to solve.

My advice is: begin to do this by yourself and come back to us if you get Maple difficulties to sove your problem.

@JoyDivisionMan 

By the way, the attached worksheet explains in more details how to get the pdf of Y = cos(X) where X  is a uniform random variable with support (-𝜋, +𝜋).
Two approaches are used: a direct construction of PDF(Y) and an indirect one by derivation of CDF(Y).

To compute the cdf of, let's say, X = sin(X) (which is obviously the same as Y) without doing all the stuff contained in the attached file (solve/allsolutions, branch identification, aggregation of partial pdf/cdf over the different branches, ...) it is much simpler to observe that sin(X) = cos(𝜋/2X) then define an auxiliary random variable A = 𝜋/2and replace everywhere it occurs  fX(something) by fA(some other thing) which is equal to fX(𝜋/2something)

restart

with(Statistics):

X := RandomVariable(Uniform(-Pi, Pi)):

Y := cos(X):

f__Y := CodeTools:-Usage( unapply(PDF(Y, y), y) )

memory used=17.78MiB, alloc change=34.00MiB, cpu time=273.00ms, real time=2.38s, gc time=0ns

 

proc (y) options operator, arrow; piecewise(y <= -1, 0, y < 1, 1/(Pi*(-y^2+1)^(1/2)), 1 <= y, 0) end proc

(1)

restart

with(Statistics):

X := RandomVariable(Uniform(-Pi, Pi)):

Y := cos(X):

F__Y := CodeTools:-Usage( unapply(CDF(Y, y), y) )

memory used=41.23MiB, alloc change=70.01MiB, cpu time=618.00ms, real time=2.02s, gc time=56.15ms

 

proc (y) options operator, arrow; piecewise(y <= -1, 0, y <= 1, (1/2)*(2*arcsin(y)+Pi)/Pi, 1 < y, 1) end proc

(2)

f__Y := unapply(diff(F__Y(y), y) , y)

proc (y) options operator, arrow; piecewise(y < -1, 0, y = -1, undefined, y < 1, 1/((-y^2+1)^(1/2)*Pi), y = 1, undefined, 1 < y, 0) end proc

(3)


FIRST WAY: BUILD DIRECTLY PDF(Y)

plot(cos(x), x=Support(X, output='range'))

 


This graph is made of two "branches" where ths cosine function has a monotone variation:
          (1)  from -Pi to 0 cosine is monotonously increasing
          (2)  from 0 to Pi  cosine is monotonously decreasing


Watchout
The inverse cosine function is the arccosine function only on intervals where cosine is monotone.
So the reciprocal of  y = cos(x) over -p..p is not x = arccos(y).
Indeed:

solve(cos(x)=y, x, allsolutions)

2*Pi*_Z3-2*arccos(y)*_B3+arccos(y)

(4)

indets((4), name) minus {y, Pi}:
about~(%):

Originally _B3, renamed _B3~:

  is assumed to be: OrProp(0,1)

Originally _Z3, renamed _Z3~:
  is assumed to be: integer

 


Choosing _Z3 = 0 gives
         

reciprocal := eval((4), _Z3=0)

-2*arccos(y)*_B3+arccos(y)

(5)


And the two branches are
         

phi[1] := unapply( eval(reciprocal, _B3=0), y);
phi[2] := unapply( eval(reciprocal, _B3=1), y);

proc (y) options operator, arrow; arccos(y) end proc

 

proc (y) options operator, arrow; -arccos(y) end proc

(6)


pdf  fX(x) of X
   

f__X := unapply( PDF(X, x), x)

proc (x) options operator, arrow; piecewise(x < -Pi, 0, x < Pi, (1/2)/Pi, 0) end proc

(7)


The pdf  fY(y) of Y is given by this formula

result := Sum('f__X'(phi[i](y)) * abs(diff(phi[i](y), y)), i=1..2)

Sum(f__X(phi[i](y))*abs(diff(phi[i](y), y)), i = 1 .. 2)

(8)

result := value(result)

result := piecewise(arccos(y) < -Pi, 0, arccos(y) < Pi, 1/(2*Pi), 0)/sqrt(abs(y^2-1))+piecewise(-arccos(y) < -Pi, 0, -arccos(y) < Pi, 1/(2*Pi), 0)/sqrt(abs(y^2-1))

(9)

f__Y := unapply(simplify(result), y):
f__Y(y);

piecewise(y < -1, 0, y = -1, 1/(2*Pi*sqrt(y^2-1)), y < 1, 1/(sqrt(-y^2+1)*Pi), y = 1, 1/(Pi*sqrt(y^2-1)), 1 < y, 0)

(10)

plot(f__Y(y), y=-1..1, title="Arcsine law")

 


ALTERNATIVE WAY: BUILD CDF(Y) AND DERIVE IT TO GET PDF(Y)


cdf  FX(x) of X
   

F__X := unapply( CDF(X, x), x)

proc (x) options operator, arrow; piecewise(x < -Pi, 0, x < Pi, (1/2)*(x+Pi)/Pi, 1) end proc

(11)


If all the fi(y) were increasing functions of y, the cdf  FY(y) of Y would write

result := Sum(chi('F__X'(phi[i](y))), i=1..2);

Sum(chi(F__X(phi[i](y))), i = 1 .. 2)

(12)


If some fi(y) are decreasing functions of y, the cdf  FY(y) of Y is defined by this formula

result := Sum(chi('F__X'(phi[i](y))), i=1..2), ``(chi('F__X'(phi[i](y))) = piecewise(diff(phi[i](y), y) > 0, 'F__X'(phi[i](y)), 1-'F__X'(phi[i](y))))

result := Sum(chi(`#msub(mi("F"),mi("X"))`(phi[i](y)))*`,`(chi(`#msub(mi("F"),mi("X"))`(phi[i](y))) = piecewise(0 < diff(phi[i](y), y), `#msub(mi("F"),mi("X"))`(phi[i](y)), 1-`#msub(mi("F"),mi("X"))`(phi[i](y)))), i = 1 .. 2)

(13)

plot([phi[1](y), phi[2](y)], y=-1..1, color=[red, blue], legend=[typeset(phi__1(y)), typeset(phi__2(y))])

 


Thus

result := ``(1 - F__X(phi[1](y))) + F__X(phi[2](y))

result*`:=`(1-piecewise(arccos(y) < -Pi, 0, arccos(y) < Pi, (1/2)*(arccos(y)+Pi)/Pi, 1))+piecewise(-arccos(y) < -Pi, 0, -arccos(y) < Pi, (1/2)*(Pi-arccos(y))/Pi, 1)

(14)

F__Y := unapply(simplify(expand~(result)), y):
F__Y(y);

piecewise(y < -1, 1, y <= 1, (Pi-arccos(y))/Pi, 1 < y, 1)

(15)


So PDF(Y) is:

diff(F__Y(y), y);

piecewise(y < -1, 0, y = -1, undefined, y < 1, 1/(sqrt(-y^2+1)*Pi), y = 1, undefined, 1 < y, 0)

(16)


Verify that fY(y) obtained using the direct way (equation (10)) is almost everywhere equal to expression (16) 

"almost everywhere" means "out of the points of null probability measure (y=±1) where (16) is undefined.
   Note that the existence of these points come from the discontinuity of PDF(X) and that they could be removed from
   formulas (10) and (16) (because of their null probability measure).

simplify(f__Y(y) - diff(F__Y(y), y));

piecewise(y = -1, undefined, y = 1, undefined, 0)

(17)

plot(
  [f__Y(y), (16)]
  , y=-1..1
  , color=[red, blue]
  , thickness=[1, 5]
  , transparency=[0, 0.7]
  , legend=[typeset('f__Y'(y)), typeset('diff(F__Y(y), y)')]
)

 


NOTE THAT MAPLE IS CAPABLE TO COMPUTE PDF(Y) DIRECTLY

PDF(cos(X), y)

piecewise(y <= -1, 0, y < 1, 1/(sqrt(-y^2+1)*Pi), 1 <= y, 0)

(18)


But I am pretty sure this is only by chance for, from what I have observed asking Maple to compute the PDF of a polynomial
function of X, Maple does not seem to identify "branches" nor to use formulas like (8) or (13).

In the cos(X) case, without no particular attention to the correctness of the formulas and by invoking wrongly a symmetry argument 
we would get:

WrongFormulaButCorrectResult := diff(simplify(2*F__X(phi[2](y))), y);

WrongFormulaButCorrectResult := diff(simplify(2*F__X(phi[1](y))), y);

WrongFormulaButCorrectResult := piecewise(y < -1, 0, y = -1, undefined, y < 1, 1/(sqrt(-y^2+1)*Pi), y = 1, undefined, 1 < y, 0)

 

piecewise(y < -1, 0, y = -1, undefined, y < 1, -1/((-y^2+1)^(1/2)*Pi), y = 1, undefined, 1 < y, 0)

(19)

V := RandomVariable(Uniform(-5*Pi/3, Pi/3)):

PDF(cos(V), y);    # Should be equal to (10) or (16)

plot(%, y=-1..1, title="Obviously wrong", titlefont=[Tahoma, bold, 14]);

piecewise(y <= 1/2, 0, y < 1, 1/(sqrt(-y^2+1)*Pi), 1 <= y, 0)

 

 
 

 

Download cosine.mw

@JoyDivisionMan 

You're welcome!

@janhardo 

I wrote "at points (a, f(a)) and (-a, -f(a)) if and only if a=0... which is impossible because f(0) is undefined" because I didn't want to write explicitely a=1 and f(a)=2... or more exactly that a was strictly larger than zero.
Indeed, if a > 0 there is no circle tangent to the curve C at points at points (a, f(a)) and (-a, -f(a)).
The only possibility for such a circle to exist is a = 0, which cannot be accepted for the function f is not defined at 0.

@acer 

Thank you for your intervention.

If this can be of any interest here is a worksheet which shows when Maple succeeds in computing the PDF and CDF of a function 𝜙 of a random variable X :

  • Maple successes are illustrated by two examples:
    • the first one is trivial because 𝜙 is linear,
    • the second one is a priori more complex because 𝜙 is not monotonous... but the "branches" have simple analytic expressions and Maple can do the job.
  • Maple failure is illustrated by the  non monotonuous function 𝜙 : X ⟼X3/3 + X 
    (note that Maple succeeds if 𝜙 : X ⟼k*X3 because  𝜙 is now monotonuous).

I know how to compute the CDF of 𝜙(X) in an ad hoc way for arbitrary functions 𝜙, but I'm not capable to know how to do it programatically in a general case.
And I guess this is the lack of such an algorithm which pushes Maple to return FAIL for almost all non monotonuous functions 𝜙 (not to speak of no-failure situations where wrong results are returned... likely for the same lack of algorithm).

When_Maple_fails_in_computing_PDF_ot_CDF.mw

@JoyDivisionMan 

I don't know if you intend to follow @acer's suggestion of submitting an SCR (Software Change Request, see the screen shot below)


But if you do it, please inform the development team that Maple 2024 and beyond must not simply reproduce the result older versions provide (2015 or instance) because this result is wrong (for your "test" case) or might be wrong in case the 
function 𝜙 : X ⟼ Y is not a one-to-one map over the whole support of X.

There is indeed an error in the algorithm which computes PDF(Y) and this is it which must be corrected.

@JoyDivisionMan

Please look my pevious answer to understand why I limit my construction to a function 𝜙 : X ⟼Y which is a one-to-one C1 map over the support of X

restart

with(Statistics):


As a rule it's almost always simpler to compute the PDF of a function of a (several) random variables by computing
firstly their CDF and differentiating the result, instead than computing directly the PDF.

Here is an example: Computing PDF(sin(X)) where X ~ Uniform(-p/2, p/2).

Note that the support of X has been chosen in order that the sine function is a one-to-one map over this support. The
case where the X Y map is not one-to-one over the support of X, for instance if X ~ Uniform(0, 2p), is a little bit more delicate to treat (see my previous answer).

X := RandomVariable(Uniform(-Pi/2, Pi/2)):

F__X := unapply(CDF(X, x), x):
F__X(x);

piecewise(x < -(1/2)*Pi, 0, x < (1/2)*Pi, ((1/2)*Pi+x)/Pi, 1)

(1)

phi := x -> sin(x);
psi := unapply(solve(phi(x)=y, x), y):

'psi'(y) = psi(y);

proc (x) options operator, arrow; sin(x) end proc

 

psi(y) = arcsin(y)

(2)

F__Y := unapply(F__X(psi(y)), y):
F__Y(y);

piecewise(arcsin(y) < -(1/2)*Pi, 0, arcsin(y) < (1/2)*Pi, (arcsin(y)+(1/2)*Pi)/Pi, 1)

(3)

f__Y := unapply(diff(F__Y(y), y), y);

proc (y) options operator, arrow; piecewise(y < -1, 0, y = -1, undefined, y < 1, 1/((-y^2+1)^(1/2)*Pi), y = 1, undefined, 1 < y, 0) end proc

(4)

plots:-display(
 Histogram(Sample(sin(X), 10^5), minbins=100, style=polygon, color="LightGray", transparency=0.5),
  plot(f__Y(y), y=-1..1, color=blue, thickness=3)
)

 
 

 

Download PDF_of_a_function_of_a_random_variable.mw

@Andiguys 

Your expressions are that complex that I am pretty sure they come from some calculus you do not give the details here, and whose results you copied-pasted into the worksheet you provide.

The reason why you get those complex expressions is likely a consequence of you doing something like
x1 := some expression containing quantities q1, .., qn
x2 := some expression containing some quantities among q1, ... ,qn plus extra quanrities r1, ..., rn
x3 := some expression containing some quantities among q1, ..., qn and r1, .., rn, plus extra quanrities s1, ..., sn
... and so on up to a quite simple expression of the form z = f(x1, x2, ..., xm)

Doing this will automatically result in an extremely complec form for z.

A better way to proceed if you want to get a concise expression for z is to write instead
z := f(x1, x2, ..., xm)
substitutions := {x1 = ?..., x2= ..., x3= ..., ...}

Now z appears in its simplest form and you don't have t to worry about how to "shorten" the expression of z.
If you run the command eval(z, substitutions), you will get an extremely complex expression that you will probably struggle to “shorten” significantly.

Here are two simple illustrations WhatYouMustNotDo.mw WhatYouMustNotDo_2.mw of what I said.

Before CAS existed, people thought twice before embarking on laborious calculations and tried to simplify the initial expressions as much as possible before undertaking these calculations, sometimes even assigning other names to sub-expressions during the calculation.
Today's CAS are so powerful that they eliminate the need to think, leading many users to want to simplify results that they have unnecessarily complicated out of laziness.
Try to avoid taking this way too.

@C_R 

I need to think about it before giving a solid answer.

Sorry.

@C_R 

Yes It is indeed a pity that Maple does not offer a function to draw (in a correct way) the cumulative distribution of discrete random variable or of a sample.

The CumulativeSumChart (see CUSUM) is essentially used in "quality control" and has nothing to do with the cumulative function (the name CumulativeSumChart can be misleading). Note that the latter is a non decreasing function while the former can have arbitrary variations. 

Cumulative Sum Charts should be used to plot stochastic process (the corresponding help page says it plots  "the cumulative sum chart for the specified data" but unfortunately gives no reference to indicate what this cumulative sum chart, which can misled those who use the CumulativeSumChart function.

For fun, look Here.mw .

@janhardo 

I vote up

@Andiguys 

For what it's worth: answer.mw

@Andiguys 

Plotting errirs are generally the simplest errors to fix because they tell you the object you want to plot is not plottable.
So simply print X_11_21 and P11 and try to determine if they are both plottable.

After correction and text location adjustments you should get this


TIP: Are both X_11_21 and P11 numeric? And if not why?   

@Andiguys 

As you mentioned the save command I guess you want to use the save/read mechanism, am I right?
So you want to save a plot structure assigned to some name in order to recover it in some other worksheet, ok?
This is completely different than generating a png/jpeg/pdf/ps file in order to include it in a paper for instance.

The  save/read  commands manipulate binary files named 'm files'.

You will find in the attached file how to use the save/read mechanism (not limited to plots because you save a name [or a sequence of names] and thus almost anything).

I nevertheless added some stuff about the use of plotsetup in case you would be interested in exporting a plot (no longer a name) into a file. I did so because your "On a new page, I want to display this image..." statement was a little bit confusing to me because save does not save an image (contrary to what plotsetup or Export do)  but the code which create this image.

Depends_on_what_you_want_to_do.mw

@acer 

The hatching was also part of your answer, even though it was not an initial request.
Aware that the were not "one of the OP's original queries" and that you had added them for your own convenience, I did not deem it necessary to send my comment to @Andiguys as well.

See you soon on another thread

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