mmcdara

3815 Reputation

17 Badges

6 years, 120 days

MaplePrimes Activity


These are replies submitted by mmcdara

@acer 

Thank you acer for this detailed answer.
I was completely unaware of such subtleties. 
There is no doubt this will help me to understand Maple a little better (even if the more I use it and vpild think I know it better and better I realize that I know it less and less)

See you soon

@acer 

Funny: to realize an elementwise operation you use map as I use to use ~.
For me, probably naively, I always thought these two ways do the same thing, at least for simple operations.

But I'm obviously wrong:

  • in the first code S is a vector and round~(S) is still of Hfloat type
  • In the second one, after conversion of S into a list, round~(S) is of integer type
    (in real situations I do not use this conversion into a list and keep operating on the vector or matiw that Sample returns... and this is why converting into integers didn't seem to me that obvious).

restart;                             

kernelopts(version);  

with(Statistics):                    

S := Sample(Binomial(10, 1/3), 3):

map(round,S)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

 

Vector[row](%id = 18446744078241542742)

(1)

round~(S);
lprint(%);

Vector[row](3, {(1) = 5., (2) = 5., (3) = 2.})

 

Vector[row](3, {1 = HFloat(5.), 2 = HFloat(5.), 3 = HFloat(2.)}, datatype = float[8], storage = rectangular, order = Fortran_order, shape = [])

 

restart;            

with(Statistics):                    

S := Sample(Binomial(10, 1/3), 3):
S := convert(S, list);

map(round,S)

[HFloat(5.0), HFloat(5.0), HFloat(2.0)]

 

[5, 5, 2]

(2)

round~(S);

[5, 5, 2]

 

integer

(3)

 

Download round.mw

restart;                             
kernelopts(version);  
with(Statistics):                    
S := Sample(Binomial(10, 1/3), 3):
map(round,S)
Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895
            Vector[row](%id = 18446744078241542742)
round~(S)
            Vector[row](%id = 18446744078241543942)
restart;            
with(Statistics):                    
S := Sample(Binomial(10, 1/3), 3): S := convert(S, list);
map(round,S)
Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895
            [HFloat(5.0), HFloat(5.0), HFloat(2.0)]
                           [5, 5, 2]
round~(S)
                           [5, 5, 2]


PS: why the insertion of a code snippet displays this ? How did you manage to do this insertion and have S displayed "in cler" ?

restart;                             
kernelopts(version);  
with(Statistics):                    
S := Sample(Binomial(10, 1/3), 3);
map(round,S)
Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895
            Vector[row](%id = 18446744078241542502)
            Vector[row](%id = 18446744078241543822)
round~(S)
            Vector[row](%id = 18446744078241544302)

@acer 

Thank you acer, you indeed raised a very interesting point.
I wasn't aware of this recision issue.

I use to use UseHardwareFloats:=false for a complete different reason (please see below) and I have observed a curious "side effect" of UseHardwareFloats.

In some algorithms one uses a Binomial trial to select the column of a matrix. Unfortunately the ouput of such trial is not an integer but an Hfloat. As converting this output into an integer is not that simple, I set UseHardwareFloats to false :

restart:
with(Statistics):
S := Sample(Binomial(10, 1/3), 1)[1];
                          HFloat(5.0)
lprint(S[1]);
HFloat(5.)[1]
round~(S);
                               5
restart:
with(Statistics):
UseHardwareFloats:=false:
S := Sample(Binomial(10, 1/3), 1)[1];
                              1.0
round~(S);
                               1

You can see that the value of S depends on the value of UseHardwareFloats, which means that either the seed used is not the same (is this seed an Hfloat?) or that the sampling algorithm doesn't proceed the same way given the value of UseHardwareFloats.

@Wes 

First question (just a complement to Carl's answer
Quantile vs solve(CDF):

restart:
with(Statistics):                   # Here I use the Statistics package
Y := RandomVariable(ChiSquare(2)):  # The correct definition of Y with this package
solve(CDF(Y, y)=0.9, y);            # a naive way
Quantile(Y, 0.9);                   # build-in procedure

                          4.605170186
                   HFloat(4.605170185988092)

Second point= What does CDF do?

restart:
with(Statistics):
Y := RandomVariable(ChiSquare(2)):  
CDF(Y, y);  # you can call CDF this way to get the Cumulative Density Function
                       /                 /  1  \\
              piecewise|y < 0, 0, 1 - exp|- - y||
                       \                 \  2  //
CDF(Y, 1);   # Personally I don't like this way which is a short of, for instance, 
             # eval(CDF(Y, y), y=1)
                                 /-1\
                          1 - exp|--|
                                 \2 /

There exits another package named Student:-Statistics, older than Statistics which does the same thing:

restart:
with(Student:-Statistics):
Y := ChiSquareRandomVariable(2):  # Here this is the correct syntax
CDF(Y, y);
CDF(Y, 1)
                       /                 /  1  \\
              piecewise|y < 0, 0, 1 - exp|- - y||
                       \                 \  2  //
                                 /-1\
                          1 - exp|--|
                                 \2 /

What I did with CDF (both packages) can be done the same way with PDF.
Also note the use of Quantile with package Student:-Statistics:

Quantile(Y, 0.9);
                   HFloat(4.605170185988092)


To sum up, there is no difference in your cas in using Statistics or Student:-Statistics.
But mixing the syntaxes could explain the problem you observed:

restart:
with(Statistics):
Y := ChiSquareRandomVariable(2):  # wrong syntax with this package
Quantile(Y, 0.9);
CDF(Y, y);
CDF(Y, 1);
                              FAIL
        piecewise(ChiSquareRandomVariable(2) <= y, 1, 0)
        piecewise(ChiSquareRandomVariable(2) <= 1, 1, 0)

Here "Y" is not a random variable for ChiSquareRandomVariable is not a function of package Statistics.
The result is thus quite logic.

The problem is that this "logic" doesn't hold for any random variable:

restart:
with(Statistics):
Y := Normal(0, 1):
Quantile(Y, 0.9);
CDF(Y, y);
CDF(Y, 1);
                   HFloat(1.2815515655447307)
                     1   1    /1    (1/2)\
                     - + - erf|- y 2     |
                     2   2    \2         /
                      1   1    /1  (1/2)\
                      - + - erf|- 2     |
                      2   2    \2       /

This is undeed quite disturbing.
Note that the names of the distributions (in both packages) can be "long names" or "short names", for instance:

restart:
with(Student:-Statistics):
Y := NormalRandomVariable(0, 1):
Quantile(Y, 0.9);
                   HFloat(1.2815515655447307)
X := Normal(0, 1):
Quantile(X, 0.9);
                   HFloat(1.2815515655447307)
U := ChiSquareRandomVariable(2):
Quantile(U, 0.9);
                   HFloat(4.605170185988092)
V := ChiSquare(2):
Quantile(U, 0.9);
                   HFloat(4.605170185988092)

and

restart:
with(Statistics):
Y := RandomVariable(Normal(0, 1)):
Quantile(Y, 0.9);
                   HFloat(1.2815515655447307)
X := Normal(0, 1):
Quantile(X, 0.9);
                   HFloat(1.2815515655447307)
U := RandomVariable(ChiSquare(2)):
Quantile(U, 0.9);
                   HFloat(4.605170185988092)
V := ChiSquare(2):
Quantile(U, 0.9);
                   HFloat(4.605170185988092)


 

@WA573 

Here is an update of my previous worksheet where I explicitely set A[0] to 6.
If you look carefully my first file you will se that I solve a system of N equations in P < N unknowns. Getting a non null solution means these equations are compatible.
Now, with A[0]=6, I have to solve a system still rectangular, but the equations are no longer compatible: trying to solve this (NxP) system gives a null answer.

If you solve all the (PxP) subsystems you can form you get different solutions.

PA_mmcdara_2.mw


If I  understand correctly your problem it seems you try to solve a differential equation.
Because dsolve doesn't  return any answer you thought to use an ansatz for u (?) 
If it is so, why don't simply use this:

eval(Eq, u=U(x));
dsolve({%, U(0)=u__0, D(U)(0)=v__0}, U(x), series, order=4);
convert(%, polynom);

Doing this explains U(x) under the form 

sum(C[m]*x^m, m=0..order-1)

 

@WA573 

Why do you say that A[0]=6? ... the instruction A[0]:=6 is commented !!!

@tomleslie  @Rakshak

What do you mean by  "eq1 and eq2 can be considered as "coupled", although "simultaneous" might be a better designation"?
Either the solution f1 of eq1 is a solution of eq2 and the solution of eq2 is also a solution of eq1, in what case eq1 anq eq2 could be said "compatible"; instead they are not and eq1 and eq2 are "uncompatible".

pdsolve([eq1, eq2]), pdsolve([eq2, eq1]), pdsolve(eq2);

both return the same solution.
But this solution is not (unless I'm mistaken) the solution of eq1: compare g1 and g2 iand the integrands o1 and o2 n the attached file. 
 

restart

eq1 := (2*(r^2+a^2*cos(theta)^2))*(M*r-(1/2)*a^2-(1/2)*r^2)*(diff(f(r, theta), r, theta))+(2*(a^2*(M-r)*cos(theta)^2-M*r^2+a^2*r))*(diff(f(r, theta), theta))

2*(r^2+a^2*cos(theta)^2)*(M*r-(1/2)*a^2-(1/2)*r^2)*(diff(diff(f(r, theta), r), theta))+2*(a^2*(M-r)*cos(theta)^2-M*r^2+a^2*r)*(diff(f(r, theta), theta))

(1)

eq2 := sin(theta)*(r^2+a^2*cos(theta)^2)*(diff(f(r, theta), theta, theta))-cos(theta)*(diff(f(r, theta), theta))*(a^2*cos(theta)^2-2*a^2-r^2)

sin(theta)*(r^2+a^2*cos(theta)^2)*(diff(diff(f(r, theta), theta), theta))-cos(theta)*(diff(f(r, theta), theta))*(a^2*cos(theta)^2-2*a^2-r^2)

(2)

eq3 := -2*(r^2+a^2*cos(theta)^2)^2*(M*r-(1/2)*a^2-(1/2)*r^2)*sin(theta)*(diff(g(r, theta), r, r))+sin(theta)*(r^2+a^2*cos(theta)^2)^2*(diff(g(r, theta), theta, theta))+(4*(-(1/4)*cos(theta)^4*a^4+a^2*r*(M-(1/2)*r)*cos(theta)^2-M*a^2*r-(1/4)*r^4))*cos(theta)*(diff(g(r, theta), theta))-2*M*sin(theta)*(diff(g(r, theta), r))*(a^2+r^2)*(cos(theta)*a-r)*(cos(theta)*a+r)

-2*(r^2+a^2*cos(theta)^2)^2*(M*r-(1/2)*a^2-(1/2)*r^2)*sin(theta)*(diff(diff(g(r, theta), r), r))+sin(theta)*(r^2+a^2*cos(theta)^2)^2*(diff(diff(g(r, theta), theta), theta))+4*(-(1/4)*cos(theta)^4*a^4+a^2*r*(M-(1/2)*r)*cos(theta)^2-M*a^2*r-(1/4)*r^4)*cos(theta)*(diff(g(r, theta), theta))-2*M*sin(theta)*(diff(g(r, theta), r))*(a^2+r^2)*(cos(theta)*a-r)*(cos(theta)*a+r)

(3)

indets(eq1, function);
indets(eq2, function);
indets(eq3, function);

{cos(theta), diff(diff(f(r, theta), r), theta), diff(f(r, theta), r), diff(f(r, theta), theta), f(r, theta)}

 

{cos(theta), diff(diff(f(r, theta), theta), theta), diff(f(r, theta), theta), f(r, theta), sin(theta)}

 

{cos(theta), diff(diff(g(r, theta), r), r), diff(diff(g(r, theta), theta), theta), diff(g(r, theta), r), diff(g(r, theta), theta), g(r, theta), sin(theta)}

(4)

f1 := pdsolve(eq1);

f(r, theta) = _F2(r)+(1/2)*(Int(_F1(theta)*(a^2*cos(2*theta)+a^2+2*r^2), theta))/(2*M*r-a^2-r^2)

(5)

f2 := pdsolve(eq2);

f(r, theta) = _F1(r)+(Int((r^2+a^2*cos(theta)^2)/((cos(theta)-1)^(1/2)*(1+cos(theta))^(1/2)), theta))*_F2(r)

(6)

pdetest(f1, eq2);
dsolve(%);
g1 := rhs(eval(f1, %));

(_F1(theta)*cos(theta)^5*a^4+(diff(_F1(theta), theta))*sin(theta)*cos(theta)^4*a^4+2*_F1(theta)*cos(theta)^3*a^2*r^2+2*(diff(_F1(theta), theta))*sin(theta)*cos(theta)^2*a^2*r^2+_F1(theta)*cos(theta)*r^4+(diff(_F1(theta), theta))*sin(theta)*r^4)/(2*M*r-a^2-r^2)

 

_F1(theta) = _C1/sin(theta)

 

_F2(r)+(1/2)*(Int(_C1*(a^2*cos(2*theta)+a^2+2*r^2)/sin(theta), theta))/(2*M*r-a^2-r^2)

(7)

pdetest(f2, eq1);
dsolve(%);
g2 := rhs(eval(f2, %));

(2*M*cos(theta)^4*(diff(_F2(r), r))*a^4*r-cos(theta)^4*(diff(_F2(r), r))*a^6-cos(theta)^4*(diff(_F2(r), r))*a^4*r^2+2*M*cos(theta)^4*_F2(r)*a^4-2*cos(theta)^4*_F2(r)*a^4*r+4*M*cos(theta)^2*(diff(_F2(r), r))*a^2*r^3-2*cos(theta)^2*(diff(_F2(r), r))*a^4*r^2-2*cos(theta)^2*(diff(_F2(r), r))*a^2*r^4+4*M*cos(theta)^2*_F2(r)*a^2*r^2-4*cos(theta)^2*_F2(r)*a^2*r^3+2*M*(diff(_F2(r), r))*r^5-(diff(_F2(r), r))*a^2*r^4-(diff(_F2(r), r))*r^6+2*M*_F2(r)*r^4-2*_F2(r)*r^5)/((cos(theta)-1)^(1/2)*(1+cos(theta))^(1/2))

 

_F2(r) = _C1/(2*M*r-a^2-r^2)

 

_F1(r)+(Int((r^2+a^2*cos(theta)^2)/((cos(theta)-1)^(1/2)*(1+cos(theta))^(1/2)), theta))*_C1/(2*M*r-a^2-r^2)

(8)

o1 := eval(op(1, numer(select(has, g1, Int))), _C1=1);
o2 := op([1, 1], numer(select(has, g2, Int)));

(a^2*cos(2*theta)+a^2+2*r^2)/sin(theta)

 

(r^2+a^2*cos(theta)^2)/((cos(theta)-1)^(1/2)*(1+cos(theta))^(1/2))

(9)

o2 := simplify(simplify(combine(o2, radical), trig)) assuming theta::real

-I*(r^2+a^2*cos(theta)^2)/abs(sin(theta))

(10)

f12 := pdsolve([eq1, eq2]):
f2, f12;

f(r, theta) = _F1(r)+(Int((r^2+a^2*cos(theta)^2)/((cos(theta)-1)^(1/2)*(1+cos(theta))^(1/2)), theta))*_F2(r), {f(r, theta) = _F1(r)+(Int((r^2+a^2*cos(theta)^2)/((cos(theta)-1)^(1/2)*(1+cos(theta))^(1/2)), theta))*_C1/(2*M*r-a^2-r^2)}

(11)

 


 

Download pde1_mmcdara.mw


3 PDEs and only 2 unknown functions (f and g)
PDEs eq1 and eq1 only depend on f and are incompatible . So either f is solution of eq1 OR (exclusive) of eq2:

pdsolve(eq1);

pdsolve(eq2);

As eq3 does not depend on your system is not a system of couple PDEs.

pdsolve(eq3);

produces no output, which means this PDE cannot be solved formally.
You can try a numerical approach but its ârameters have to be geven numerical values, AND you must provide boundary condifions.
 

@MaPal93 

I never use the Document Mode with 2D inputs, never, never!
IMO it is too sensitive to typing errors which can insert hidden characters (browse que Question section here to make your own idea about the pros and cons of this inpit mode).

Here is your file with classical Maple input inserted (maube less beautiful but it works correctly)
sub_matrices_toplot_mmcdara.mw

I know I'm not answering your question, but I hate this input mode so much that I don't want to waste time trying to figure out where your error comes from.

@Carl Love 

My answer was based on my understanding, of what I thought the OP wanted.
Regarding the dualaxisplot, I totally agree with you: in this case, it looks more like comparing two graphs and a simple display seems more appropriate.

@MaPal93 

I have thought it would be obvious, sorry:

plots:-display(`<|>`(TwoMat(A_jk, A_ki), TwoMat(A_jk, A_ki, scale = true)));

# Or if you just want unscaled matrices:
TwoMat(A_jk, A_ki)

Maybe you are puzzled by the plot if A_ki ???
In this case look to this matrix or do 

plots:-matrixplot(A_ki, heights = histogram)

to understand what happens

@dharr 

  1. I need to think again to that.
    "minimizing only finds one of the solution": right, but I believe that there is (maybe I'm wrong) a single set of values for the lambda(s) and an infinity for the mu(s); something as if the solution was a variety of dimension 1 in a space of dimension 8 (???).
    Observation: starting from different initial points didn't give me different results for the lambda(s), only for the mu(s).
     
  2. Observation: It's clear that 6 equations depend only on the lambda(s), but applyng fsolve on this sub system gives no solution whatever the value of Digits or the initial point (for those I used).
     
  3. By the way V=~0 is not a solution:
V=~0;
eval(Eqp,%);

Thanks for all your comments

@MaPal93 

Optimization:-NLPSolve is faster and more stable than fsolve:

matrix_inverse_mmcdara.mw

Splitting the system Eqp in two sub-systems reveals that its solution verifies 

{mu__jk = -mu__ki, mu__ki = mu__ki}

 

The dyadic product can be realized in a simple way without using LinearAlgebra:-KroneckerProduct:

N := 3: 
A := Vector(N, symbol=a): 
B := Vector(N, symbol=b): 
M1 := LinearAlgebra:-KroneckerProduct(A, B^+): 

# alternative way

M2 := Matrix(N$2, (i, j) -> A[i]*B[j]):

Which method is the fastest?

 

restart:

N := 3:
A := Vector(N, symbol=a):
B := Vector(N, symbol=b):
LinearAlgebra:-KroneckerProduct(A, B^+):
M := Matrix(N$2, (i, j) -> A[i]*B[j]):

%%, %

Matrix(3, 3, {(1, 1) = a[1]*b[1], (1, 2) = a[1]*b[2], (1, 3) = a[1]*b[3], (2, 1) = a[2]*b[1], (2, 2) = a[2]*b[2], (2, 3) = a[2]*b[3], (3, 1) = a[3]*b[1], (3, 2) = a[3]*b[2], (3, 3) = a[3]*b[3]}), Matrix(3, 3, {(1, 1) = a[1]*b[1], (1, 2) = a[1]*b[2], (1, 3) = a[1]*b[3], (2, 1) = a[2]*b[1], (2, 2) = a[2]*b[2], (2, 3) = a[2]*b[3], (3, 1) = a[3]*b[1], (3, 2) = a[3]*b[2], (3, 3) = a[3]*b[3]})

(1)

N := 1000:
A := Vector(N, symbol=a):
B := Vector(N, symbol=b):
CodeTools:-Usage( LinearAlgebra:-KroneckerProduct(A, B^+) ):
CodeTools:-Usage( Matrix(N$2, (i, j) -> A[i]*B[j]) ):

memory used=91.63MiB, alloc change=339.64MiB, cpu time=3.36s, real time=2.85s, gc time=1.65s

memory used=45.79MiB, alloc change=7.63MiB, cpu time=876.00ms, real time=615.00ms, gc time=356.94ms

 

 

Download DyadicProduct.mw

As it sometimes happens, the built-in procedure is neither the fastest nor the most parsimonious in terms of memory allocation.

 

@Dkunb 

A faster algorithm to plot the level curve f(a, b)=C.
(same toy problem than in my previous answer).
 

restart:

with(plots):

hA := 0.2:
hB := 0.2:

alpha := Array([seq(0..10, hA)]):
beta  := Array([seq(0..10, hB)]);
PA := numelems(alpha);
PB := numelems(beta);

beta := Vector(4, {(1) = ` 1 .. 51 `*Array, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

51

 

51

(1)

NN := Matrix(PA, PB, (i, j) -> cos(alpha[i]/3+beta[j]/3)*sin(alpha[i]*beta[j]/9))

NN := Vector(4, {(1) = ` 51 x 51 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(2)

phi := (a, b) -> [(1-a)*(1-b), (1-a)*b, a*b, a*(1-b)]

proc (a, b) options operator, arrow; [(1-a)*(1-b), (1-a)*b, a*b, a*(1-b)] end proc

(3)

C  := 0.2:
LC := NULL:
for i from 1 to PA-1 do
  for j from 1 to PB-1 do
    nodes  := NN[i, j], NN[i, j+1], NN[i+1, j+1], NN[i+1, j]:
    bounds := (min..max)(nodes):
    if verify(C, bounds, interval) then
       loc := add(phi((a-alpha[i])/hA, (b-beta[j])/hB) *~ [nodes]):
       s   := solve(loc, b):
       LC  := LC,
              contourplot(
                add(phi((a-alpha[i])/hA, (b-beta[j])/hB) *~ [nodes]),
                a=alpha[i]..alpha[i+1],
                b=beta[j]..beta[j+1],
                contours=[C]
              )
    end if:
  end do:
end do:

display(LC, view=[0..10, 0..10])

 

 


 

Download LevelCurve_2D.mw

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