mmcdara

2000 Reputation

14 Badges

4 years, 321 days

MaplePrimes Activity


These are answers submitted by mmcdara

Here is a way:

restart:
with(plots):
rhsz := -y(t)*abs(y(t)):
sys2 := {
  diff(z(t), t) = rhsz
  ,diff(y(t), t) = z(t)
  ,y(0) = 1.0, z(0) = 0.
  ,d3(t)=diff(rhsz, t)  
  ,d3(0)=__d3             
  ,d4(t)=diff(d3(t), t)
  ,d4(0)=__d4
}:

tmax := 10:
dsol2 := dsolve(sys2, numeric, parameters=[__d3, __d4], range = 0 .. tmax); 

dsol2(parameters=[0, 0]):  # the values of the parameters to use is an issue
dsol2(tmax);


HigherOrderDerivatives.mw

 

I don't know what is to be qualified as a bug in Maple.
Sometimes things that I consider myself a bug are considered not to be a bug by others and reciprocally.
When things go wrong with piecewise (which I don't necessarily consider a bug but a strategy that Maple doesn't use), I usually convert the piecewise component to Heaviside and back to the piecewise form. 

f := t -> piecewise(0 <= t, 1, t < 0, t - 1);
g := unapply(convert(f(t), Heaviside), t);
G := x -> int(g(t), t = 1 .. x):
G(x);
convert(%, piecewise)

         /            1  2                                    \
piecewise|x < 0, -1 + - x  - x, x = 0, undefined, 0 < x, x - 1|
         \            2                                       /

 

Many things to fix:

  1. declare delta as local,
  2. as delta is a free parameter, use the option parameters=[delta] in dsolve
  3. instanciate res1 (=assign a value to delta) before plotting
  4. use plots:-odeplot instead of plot or plot3d
restart:
local delta:
phi:=0:lambda:=0.1:N:=5:M:=sqrt(N*(N+1))*exp(I*phi):omegap:=10:
var:={n(t),u(t),z(t)}:
dsys:={diff(n(t),t)=-2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda)+((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),diff(u(t),t)=-2*(1-I*delta)*u(t)+2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,diff((z(t),t))=-2*(1+I*delta)*z(t)+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*conjugate(M)};
 

pars := convert(indets(dsys, name) minus {t}, list);
                            [delta]

res1:=dsolve(dsys union {n(0)=0,u(0)=0,z(0)=0},numeric, parameters=pars):

# exemple
res1(parameters=[1]):
TMAX := 1:
res1(TMAX);

# Exemple 1

plots:-odeplot(
  res1,
  [t, n(t)],
  t=0..1
)

A few example (without cosmetics) are given in the attached file (with a liitle add concerning the option range)
Dsys.mw

A suggestion : do not assign values to lambda, N and omegap but include them as parameters (the  command pars := ... will do the job automatically)

IGNORE THIS ANSWER, EVERYTHING IN IT IS WRONG
mmcdara


In your set of solutions several consecutive ones are wrong, for instance   4.155071365 and 4.158473414 (one of them is a spurious one).

This issue, and the fact some solutions are skipped, comes from the fact that f(x) is always positive and that zeroes are points where f(x) tangents the x axis and not points where f(x) cuts it.
Here is a way to overcome this problem (still using NextZero).

I use an auxiliary function f__eps(x) equal to f(x)-eps to predict the location a the next zero of f(x), and then start from this predictor to find the zero of f(x).
Obviously the zeros of  f__eps(x) are the points where the graph cuts the x-axis.

restart: 
Digits := 10;
f      := proc (x) options operator, arrow; x*(1+1.0001*sin(x^2)) end proc:
eps    := 1:
f__eps := proc (x) options operator, arrow; x*(1+1.0001*sin(x^2))-eps end proc:


sol := NULL:
s := 1.0:
s := RootFinding:-NextZero(f, s): 
sol := sol, s:
printf("%3d  %2.8f\n", 1, s):

for j from 2 to 100 do 
  s__eps := RootFinding:-NextZero(f__eps, s): 
  s      := RootFinding:-NextZero(f, s__eps): 
  sol    := sol, s:
  printf("%3d  %2.8f\n", j, s):
end do:

Here is the result for the 100 first zeros : none of them is skipped and no spurious pairs are found (the blue lines locate the zeros).


Predictor_Corrector.mw  (updated)

 

One can compute HazardRate in two ways:

  1. In this simple case it's immediate to obtain that Z =1/2*X1 + 1/2*X2 is Triangular(0, 1/2, 1)
    Z := RandomVariable(Triangular(0, 1, 1/2)):
    
    HazardRate(Z, 0.4);
    HazardRate(Z, 2/5);
    
                              2.352941176
                                   40
                                   --
                                   17
    

     

  2. But this is cheating in some sense.
    The problem comes from the fact that setting Y :=1/2*X1 + 1/2*X2  does not give Y the attributes of a RandomVariable.
    This is simply shown in the attached file:
    att := attributes(X1);
      protected, RandomVariable, _ProbabilityDistribution
    exports(att[3]);
    Conditions, ParentName, Parameters, CDF, HodgesLehmann,  InverseSurvivalFunction, Mean, Median, MGF, Mode, PDF, Quantile, RousseeuwCrouxSn, Support, Variance, RandomSample, RandomSampleSetup, RandomVariate, MaximumLikelihoodEstimate
    
    att := attributes(Y);
    exports(att[3]);
    Error, invalid subscript selector
    

The idea (no cheating here) is to define a  RANDOM VARIABLE Ytrue whose PDF and CDF are those derived from the definition of Y:
 

Y     := 1/2*X1 + 1/2*X2:
dist  := Distribution(PDF = unapply(PDF(Y, t), t), CDF = unapply(CDF(Y, t), t)):
Ytrue := RandomVariable(dist):
HazardRate(Ytrue, 0.4);
HazardRate(Ytrue, 2/5);
                       
                          2.352941176
                               40
                               --
                               17

Note that attributes(Ytrue) strangely keeps retuning an error ????

You can also define Ytrue in a more detailed way, for instance 

dist  := Distribution(
            PDF = unapply(PDF(Y, t), t), 
            CDF = unapply(CDF(Y, t), t), 
            HazardRate = unapply(PDF(Y, t)/(1-CDF(Y, t)), t)
         );
Ytrue := RandomVariable(dist):

# which will enable you to get a formal expression for the hazard rate

HazardRate(Ytrue, y);


HazardRate.mw

To return to your answer to vv, even if Statistics is a good package (I agree with vv), it still suffers from several flaws (I have already pointed out some of them in the past).

And... YES, I consider you put your finger on a bug, on of those serious bug where Maple returns an answer instead of modestly saying "I can't compute that"

in the plot before the error :

  • copy the part of the command < ...>
  • type CTRL+J
  • past in the new block, you get this
    `<,>`(seq*(eval(mysolh, [z = res[i, 1], h1(z) = res[i, 2], h2(z) = res[i, 3], h3(z) = res[i, 4]]), i = 2 .. op([1, 1], res)))

    This shows a spurious character (*)  after seq.

  • Remove the lnk character between seq and (
    Again copy-paste to verify you get this

    `<,>`(seq(eval(mysolh, [z = res[i, 1], h1(z) = res[i, 2], h2(z) = res[i, 3], h3(z) = res[i, 4]]), i = 2 .. op([1, 1], res)))

     

Now the plot command is almost correct because  res[2..-1, 1] is a real valued vector while <...> is a complex valued vector.
It's now up to you  to fix this (maybe by writting Re(<...>) ?)

Simple_Plot_correction.mw

PS : try not to use the 2D input mode to avoid such kind of problems


using MathML:
 

ih := `#mi(h))`:
plot(ln(h), h=1..10, labels=[h, ln(ih)])

... IMA the more elegant way to answer the question
 

restart
with(geom3d):

point(__O, 0, 0, 0);
point(A, 6, 0, 2);
point(B, 0, 3, 2);
point(C, 0, -5*cos(Pi/3), 5*sin(Pi/3));

# Angle BAC

line(L__AB,[A, B]);
line(L__AC,[A, C]);

FindAngle(L__AB, L__AC);
`#mover(mo("BAC"),mo("&#x5e;"))` = evalf(%*180/Pi):

# Angle OAC

line(L__AO,[A,__O]);
line(L__AC,[A, C]);

FindAngle(L__AO, L__AC);
`#mover(mo("OAC"),mo("&#x5e;"))` = evalf(%*180/Pi)

angles.mw

Simply

fourier(t^n*exp(-2*t)*Heaviside(t), t, f) assuming n > 0:
G := unapply(%, n):     
plot([seq(abs(G(n)), n=0..4)], f=0..10, legend=''G''~([$0..4]))

 

 The key is that CharacteristicPolynomial always returns the eigenvalues in the same order

restart:
with(LinearAlgebra):

lambda := < [solve(CharacteristicPolynomial(Sy, t))] >:
`<|>`(seq(LinearSolve(Sy-DiagonalMatrix([L$3]), <0,0,0>), L in lambda)):
lambda, eval(%, [indets(%)[]]=~[1$3]);

 

 

Here is an answer to your first question only.
It uses two classical results about Legengre polynomials

# Legendre's differential equation

LDE := (1-x^2)*Diff(P(n, x), x$2) = 2*x*Diff(P(n, x), x) - n*(n-1)*P(n, x);

# 1st Derivative (classical result)

FD := Diff(P(n, x), x) = n/(1-x^2) * (P(n-1, x) - x*P(n, x));

# Expression of Diff(P(n, x), x, x):

map(collect, isolate(eval(LDE, FD), Diff(P(n, x), x$2)), P(n, x))

When P(n, x) is written P(n, cos(theta)), the final result is of the form

Diff(P(n, cos(theta)), theta$2) = a(n, theta)*P(n-1, cos(theta)) + b(n, theta)*P(n, cos(theta))

Where a(n, theta) and b(n, theta) are given in closed forms in the attached file

legendre.mw

PS: hope I wasn't mistaken...

A workaround

restart:
alias(u=u(x,y)):
int(diff(u,x)^2+u*diff(u, x$2), x):  # doesn't give the expected result

# workaround
with(IntegrationTools):
J := Expand(Int(diff(u*diff(u, x), x), x));  # note the use of Int instead of int
K := Parts(op(2, J), u);
sol := op(1, J) + K

integration.mw

Seems simpler

map(k -> v[.., k]/norm(v[.., k], 2), [$1..3])

 

In four steps:

  1. Do as Kitonum said (export as image.jpg for instance)
  2. Use some digitalizer (on Mac I use Engauge Digitizer, see here http://digitizer.sourceforge.net ) to gigitalize the curve (quite easy, the main problem is a correct definition of the axes.
    It's quite fast (less than 2 minutes for the result shown below).
    Export the digitalization as digitalized.csv
  3. Open digitalized.csv and replace the decimal separator (comma) by a dot
    Save the result (same name for instance)
  4. Read digitalized.csv in your worksheet (you can use CurveFitting:-Spline to create a finer digitalization):
    restart:
    file := cat(currentdir(), "/Desktop/digitalized.csv"):
    M := ImportMatrix(file):
    plots:-pointplot(M[2..-1], style=line, size=[800, 400])
    

    Excerpt from M

    numelems(M[..,1]):
    convert(M[100..102,..], listlist)
    
    [[202.3029, 204.884], [203.3591, 205.464], [203.6805, 205.47]]
    
    


    Result:


     

This can be done using MathML coding/
Example

restart

 # example 1

`#msup(mo(u),mo(2))` := 3;
`#msub(msup(mo(u),mo(2)),mo(h))` := 4;
3*`#msup(mo(u),mo(2))` + a*`#msub(msup(mo(u),mo(2)),mo(h))`;

# example 2
 
`#msup(mo(u),mo(2))` := '`#msup(mo(u),mo(2))`':
`#msub(msup(mo(u),mo(2)),mo(h))` := ' `#msub(msup(mo(u),mo(2)),mo(h))`':

f := (x, y) -> 1/log(x)+y^3:
f(`#msup(mo(u),mo(2))` , `#msub(msup(mo(u),mo(2)),mo(h))`);

But manipulating such names is not that easy.
Using aliases can help simplifying all this:

restart

# alias definitions

alias(`#msup(mo(u),mo(2))`=__x, `#msub(msup(mo(u),mo(2)),mo(h))`=__y);
       #msup(mo(u),mo(2)), #msub(msup(mo(u),mo(2)),mo(h))

# example 1

3*__x+a*__y;

# example 2
​​​​​​​
f := (x, y) -> 1/log(x)+y^3:

f(__x, __y);  # smart output
lprint(%);    # detailed output

 

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