mmcdara

4019 Reputation

17 Badges

6 years, 177 days

MaplePrimes Activity


These are replies submitted by mmcdara

 @zenterix   

When is (-2)^x real ?
The simplest way (IMO) to see this is to rewrite (-2)^X in trigonometric form:

convert(convert((-2)^~X, exp), trig);

   (cosh(ln(2) X) + sinh(ln(2) X)) (cos(Pi X) + I sin(Pi X))

# ==> (-2)^x  is real iif sin(X*Pi) = 0:

solve(sin(X*Pi) = 0, allsolutions);
               _Z1

about(%);
Originally _Z1, renamed _Z1~:
  is assumed to be: integer

# thus (-2)^x is real iif is an integer

The trigonometrixc representation also shows that (-2)^X is imaginary iif X=Z+1/2 where Z is an integer.
If X <> (Z::integer)/2 then (-2)^X is complex (non zero real and imaginary parts).

@Carl Love Not that surprising once (-2)^X is written in trigonometric form: the expression

(cosh(ln(2) X) + sinh(ln(2) X)) (cos(Pi X) + I sin(Pi X))

represents the equation of a spiral in the complex plane. Thus the smoothness of the curve.

@Carl Love 

You guess well, I only search trace within the NLPSolve and Optimization help pages.
I probably thought that a tracing command that would support many different procedures might be a utopia.I am going to quickly look at how to use it and what I can do with.
 

@x Im tc 

The formula (10) is correct as soon as Sa > Sd > 0
Maybe this is no longer the case if Sa > Sd (I did not check this point).
The question is: What do you mean by "generic" ?

@WA573 

As

2.*sin(phi-1.550000000)

takes values in the range -2..2, then f defined by

f := (phi, h) -> sqrt(2.*h+2.*sin(phi-1.550000000))

takes real values iif 2*h > 1.

The plots you get are misleading: they are not the plots of the 6 (4  more precisely) functions, but those of their real parts.
Do this to convince you:

plot(Re~([s]), phi=-2*Pi..2*Pi, size=[800,500], title="Real part");
plot(Im~([s]), phi=-2*Pi..2*Pi, size=[800,500], title="Imaginary part");

The first one is a better representations of yours (look as in your plots the curvves are not entirely drawn close to the phi axis, and hox they are now).
The second plot proves the imaginary part is not null.

So it's up to you to know what you want:

  • The functions all real valued?
  • Complex functions and thus real and imaginaty plots (or other types of plots, look the plots help page).
    Here is an example:
     

    restart;

    Take m = 1 and renane "(&PartialD;phi)/(&PartialD; xi)" to y:

    1/2*y^2 - (sin(phi-1.55)) = h:
    solve(%, y):
    f := unapply(%[1], phi, h):

    s := seq(f(phi,h), h in [-1,0,1]), seq(-f(phi,h), h in [-1,0,1]):


    col := () -> ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]):

    plots:-display(
      seq(
        plot([(Re, Im)~(s[k])], phi=-2*Pi..2*Pi, color=[col()$2], linestyle=[1, 3],legend=[typeset('Re'(k)), typeset('Im'(k))]),
        k=1..6
      ),
      size=[600,500]
    );

     

     


  • Download mww_mmcdara.mw

@x Im tc

Of course such formulae exist.
Have you read my answer (formula 10) ?
If so, why are you asking @tomleslie? Is there something wrong with what I did? Did I not already answer your question?


@tomleslie

You doubt a generic formula exists?
But it does, just read carefully the answer(formula 10) I posted 4 hours before yours, instead of simply using the same couple (Sa, Sd) and deliver a pale copy. 

@tomleslie 
Good point, I forgot to mention that I was using Windows.
I vote up

@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.

2 3 4 5 6 7 8 Last Page 4 of 91