2000 Reputation

14 Badges

4 years, 321 days

MaplePrimes Activity

These are answers submitted by mmcdara

Here is a way:

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)  
  ,d4(t)=diff(d3(t), t)

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


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):
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
local delta:

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

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

# exemple
TMAX := 1:

# Exemple 1

  [t, n(t)],

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

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


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.

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).  (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. 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
    Conditions, ParentName, Parameters, CDF, HodgesLehmann,  InverseSurvivalFunction, Mean, Median, MGF, Mode, PDF, Quantile, RousseeuwCrouxSn, Support, Variance, RandomSample, RandomSampleSetup, RandomVariate, MaximumLikelihoodEstimate
    att := attributes(Y);
    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);

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);

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(<...>) ?)

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


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__AC,[A, C]);

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


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


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

PS: hope I wasn't mistaken...

A workaround

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

# workaround
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

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 ) 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):
    file := cat(currentdir(), "/Desktop/digitalized.csv"):
    M := ImportMatrix(file):
    plots:-pointplot(M[2..-1], style=line, size=[800, 400])

    Excerpt from M

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



This can be done using MathML coding/


 # 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:


# 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


# 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