Deliberately creating noisy data is not uncommon.

Generate points along the sin curve, stored in a Vector of x values and a Vector of y values. Generate a noise signal from a random distribution, generating a sample of noise values in a Vector. Add together the y-values Vector and the noise Vector. You are done. View with ScatterPlot, as one way to handle noisy experimental data.

So, just use Statistics:-Sample to generate noise, and add noise to pure signal. You can code it, I'm sure.

acer

The following has more visible brackets (which you didn't ask for), but the final result of 16 is accesible from `h` after the assignment, and it doesn't require that you encode the formula more than once when typing it in.

restart:

twostep:=proc(expr::uneval)
   (convert(expr,name) = 
    subsindets(expr,And(name,satisfies(t->type(eval(t),constant))),
               z->convert(eval(z),name)))
   = eval(expr):
end proc:

a:=1:
b:=3:
c:=5:

h:=twostep(a+b*c);

                     (a+b*c = 1 + 3 5) = 16

rhs(h);

                               16

lhs(h);

                        a+b*c = 1 + 3 5

h:=twostep(a+b.c);

                   (a+b . c = 1 + 3 . 5) = 16

twostep( sin(exp(a*x)+b)/c );

  /                    sin(exp(1 x) + 3)\   1                
  |sin(exp(a*x)+b)/c = -----------------| = - sin(exp(x) + 3)
  \                            5        /   5                


restart:

twostep:=proc(expr::uneval)
   (convert(expr,name) = 
    subsindets(expr,And(name,satisfies(t->type(eval(t),constant))),
               z->``(eval(z))))
   = eval(expr):
end proc:

a:=1:
b:=3:
c:=5:

h:=twostep(a+b*c);

                  (a+b*c = (1) + (3) (5)) = 16

rhs(h);
                               16

lhs(h);
                     a+b*c = (1) + (3) (5)

twostep( sin(exp(a*x)+b)/c );

/                    sin(exp((1) x) + (3))\   1                
|sin(exp(a*x)+b)/c = ---------------------| = - sin(exp(x) + 3)
\                             (5)         /   5                

See also here, and maybe also this or that craziness.

acer

You're asking for an absolute conversion, and getting a relative conversion.

In relative terms, if you increase a quantity by a single unit Celcius then the corresponding increase in terms of units Kelvin should not be done on the absolute scale! Ie, it is not a bug. Eg. 55*Unit(K)+Unit(degC) should only add one Unit(K) to get the result.

The `convert` command offers you the choice. You can use convert,units for a relative scale conversion, or you can use convert,temperature for an absolute scale conversion (which you seem to be wanting).

I posit that using an absolute scale conversion for arithmetic is not usual.

You can replace your 11*Unit(degC) with the call convert(11*Unit(degC), temperature, K) in your sequence of names=values.

Of course, you could alternatively adjust literally by -273.15, but for other units like degF, etc, you might prefer it error free and programatic. Speaking of avoiding human error, your other "correct" worksheet adjusted by 273, but it ought to be 273.15.

question to experts: Would it make sense for there to be a Units `context` of `temperature`, so that one could enter something like Unit(degC[temperature]) and have conversions be done on an absolute scale after issuing UsingContexts('temperature')? It could get messy, because one wouldn't want 1*Unit(degC) + 1*Unit(degC) to come out wrongly as 548.3*Unit(K). Conversion would have to be done only after combining. Is it unworkable, on account of some other problematic cases? Would it make any sense for Units:-AddUnit to allow a procedure instead of an algebraic expression for the 'conversion'?

acer

You might find this example of interest. And you can adjust it, so that the coloring scheme is based on the orientation as opposed to simply the x-y-z position. You could then make an animation from this, giving the effect of coloring that changes with the orientation.

Getting it to work using the built-in 3D plot rotation mechanism, is trickier. During other recent investigations of my own, I've been lamenting the lack of programatic access to the current orientation of a manually rotated plot. I would like to see the orientation angles be (new) properties of a Plot Component, that could be used with the GetProperty and SetProperty commands of the DocumentTools package.

acer

One does not need to assign, in order to use the value found.

The result from Optimization:-Minimize is a list of two things, and either can be accessed by simply indexing into that list (ie. 1st, or 2nd element).

sol := Optimization:-Minimize( (x-5)^2-x, x=-10..10 );
              sol := [-5.2500000000000, [x = 5.50000000000000]]

sol[2];
                           [x = 5.50000000000000]

sol[1];
                              -5.2500000000000

eval( (x-5)^2-x, sol[2] );
                              -5.25000000000000

eval( x, sol[2]);
                              5.50000000000000

acer

The initial parsing by HTTP:-Get in Maple 12 seems to be stripping out the `?`, retaining it for neither (what it considers as the table entry indexed by) "urlpath" nor "query". That's what it looks like, to me, in the Maple debugger.

So I tried repeating the "?" mark, wondering whether the url parser would not notice and so only strip one instance of the symbol and retain the duplicate. It seems to work, for this example. That is, in Maple 12.02 (Windows 7, 64bit),

HTTP[Get]("www.mapleprimes.com/recent/all??page=2");

acer

restart:

Digits:=15:

b:=-95: `ε`:=1/4: h:=.69: n:=9.2: om:=.7:

ode1 := diff(Omega(z),z)
        +(Omega(z)*(1-Omega(z))*(3-2*sqrt(Omega(z))*(1+z)/n)
        +Omega(z)*(1-Omega(z))*`ε`
        *b*(1+z)*(diff(phi(z),z)))/(1+z)=0:
ode2 := diff(phi(z),z)+sqrt((2*(1+z))*sqrt(Omega(z))/(3*n)
             +(1-Omega(z))*(1+z)*b*`ε`
             *(diff(phi(z),z))/(3*Omega(z)))/((1+z).H(z))=0:
a := solve(ode2, diff(phi(z),z)):
ode22 := diff(phi(z),z) = a[1]:
ode3 := diff(H(z),z)+H(z)*((3/2)*Omega(z)
        -(1+z)*Omega(z)^(3/2)/n+(1/2*(1-Omega(z)))*(1+z)*b*`ε`
        *(diff(phi(z),z))-3/2)/(1+z)=0:

sys := {ode1,ode22,ode3}:
ics := {H(0)=.69, phi(0)=0, Omega(0)=.7}:

sol := dsolve(`union`(sys, ics), numeric, relerr=1e-5,
              output=listprocedure):
HH := subs(sol, H(z)):

st := time():
evalf(Int(z->1/HH(z), 0..1091.3, epsilon=1e-5)):
evalf[7](%);
                            4.566725

time()-st;
                             2.481

sol := dsolve(`union`(sys, ics), numeric, relerr=1e-5,
              output=listprocedure, stiff=true):
HH := subs(sol, H(z)):

st := time():
evalf(Int(z->1/HH(z), 0..1091.3, epsilon=1e-5)):
evalf[7](%);
                            4.566721

time()-st;
                             9.906

acer

That's what the fsolve command will do with it.

> restart:

> infolevel[fsolve]:=6:

> fsolve({4*y^2+4*y+52*x=19,169*x^2+3*y^2+111*x-108*y=10});

fsolve/sysnewton: trying multivariate Newton iteration
fsolve/sysnewton:
guess vector Vector(2, [6.5647247568434774,-9.1035796650957585])
fsolve/sysnewton: norm of errors: 9851.1068410628682
fsolve/sysnewton: new norm: 2375.9292767215670
fsolve/sysnewton: iter = 1 |incr| = 9.8322 new values x = 3.0420 y = -2.7942
fsolve/sysnewton: new norm: 499.41121284235334
fsolve/sysnewton: iter = 2 |incr| = 6.0295 new values x = 1.5984 y = 1.7918
fsolve/sysnewton: new norm: 172.63168891486316
fsolve/sysnewton: iter = 3 |incr| = 2.8820 new values x = .66908 y = -.16092
fsolve/sysnewton: new norm: 18.522740642845857
fsolve/sysnewton: iter = 4 |incr| = .77961 new values x = .35166 y = .30127
fsolve/sysnewton: new norm: .21835042965281877
fsolve/sysnewton: iter = 5 |incr| = .13242 new values x = .32250 y = .40453
fsolve/sysnewton: new norm: 0.112351888740128e-3
fsolve/sysnewton: iter = 6 |incr| = 0.85011e-3 new values x = .32168 y = .40450
fsolve/sysnewton: new norm: 0.6995170e-11
fsolve/sysnewton: iter = 7 |incr| = 0.93993e-6 new values x = .32168 y = .40450
fsolve/sysnewton: new norm: 0.917e-15
fsolve/sysnewton: iter = 8 |incr| = 0.98079e-13 new values x = .32168 y = .40450

              {x = 0.3216832637, y = 0.4044985193}

> solve({4*y^2+4*y+52*x=19,169*x^2+3*y^2+111*x-108*y=10},Explicit,AllSolutions):

> evalf([%]);

            [{x = -2.054245812, y = 5.130736680}, 

              {x = 0.3216832636, y = 0.404498520}, {

              x = 0.4402457709 + 1.697384362 I, 

              y = -3.767617600 + 3.376465579 I}, {

              x = 0.4402457709 - 1.697384362 I, 

              y = -3.767617600 - 3.376465579 I}]

I find your second example to be ambiguous. Alongside your use of implicit multiplication, I am not sure what you mean by "e^-1x". Do you mean exp(-1)*x, or exp(-x)? But you could call fsolve on the second system in a similar way to the first example.

I wasn't sure whether you wanted all the real numeric solutions, or actually needed to use Newton's method specifically to find a (any) root. There are other way to try and compute all roots, say within a given finite range.

acer

For the commandline interface running on a color-capable Linux/Unix terminal, the following prints as green whatever substring is first processed by the procedure `makegreen`. Following the green, is the code to switch to white (since that is likely the default color, or what one started with). Adjust accordingly.

makegreen:=s->sprintf("\e[32m%s\e[00m",s):

printf("%s\n", makegreen("Is this green?"));

printf("Good %ld morning. %s %s\n",
       42, makegreen("Is this green?"), "And is this not?");

And so on.

It is using ANSI escape codes, so it won't work in any shell that doesn't support interpretation of ANSI escape codes. I don't know what support for that there is in cygwin, or the old DOS shell running (or being faked) in modern MS-Windows, or in Windows Powershell.

It probably works on a OSX shell that supports ANSI color, too.

It looks like `fdiscont` has more success than symbolic `discont` here.

Try these,

plot(F(x),x=-1..27,y=0..0.5,gridlines=true,discont=true,usefdiscont);

plot(F(x),x=-1..27,y=0..0.5,gridlines=true,discont=[usefdiscont]);

The following also gets it right. I suspect that's because `discont` needs an expression (like your unevaluated F(x) as the first argument to `plot`). So, since `discont` doesn't accept an operator, like just F instead of F(x), then plot has to fall back to using `fdiscont`.

plot(F,-1..27,0..0.5,gridlines=true,discont=true);

The first pair were about 2-3 times faster than the last one.

It might look more legible (to some people), but creating an operator like `F` only so that it can ever get used as a function call like F(x) is something that often leads to confusion because it muddles unnecessarily with the distinction between operator and expression. You didn't encounter it here (maybe because the `Statistics` author went gone to extra pains to anticipate it) but this confusion often results in an error message, when the operator is prematurely evaluated at symbolic `x` while it is only prepared to handle numeric values for parameter x.

Another possibility that works is,

plot(CDF(X,x),x=-1..27,y=0..0.5,gridlines=true,discont=[usefdiscont]);

where once again, the piecewise symbolic result from CDF(X,x) may not get best result from symbolic `discont`.

acer

Wrap the other `plot` calls (or plots:-display calls) each in their own `print` function call, within the procedure.

For example,

proc(..)
  local P, Q;
  .
  .
  print( plot(...) ); # amd/or other  lines similar
  .
  .
  P := plot(...);
  Q := plot(...); # etc, or in a loop
  .
  .
  print( P );
  print( Q ); # etc
  .
  .
  return somethingelse;
end proc:

acer

If the GUI's displayprecision isn't supprted within the legend of an inlined Standard GUI plot then you could still control it manually using either `evalf` or possibly some combination of printf (or maybe even sprintf and convert/name).

Eg,

      'legend' = [ seq( 'D(t)' = evalf[5](legendList[k]) , k = [1,2] ) ]

Note the the GUI's displayprecision is pretty much entirely a cosmetic effect on its display. Since it doesn't affect the underlying value then, even if it worked within a legend, the export drivers would likely still see the original value of longer digits. So actually manipulating the digits of the value, using evalf or s/printf, could give more joy.

acer

Using 64bit Maple 15.01 on Windows 7,

restart:

plotsetup('ps', 'plotoutput'="c:/temp/P.ps",
          'plotoptions'=`noborder,landscape,height=300,width=700`):

P:=plot(thing^2, thing=0..1, labels=[thing,stuff],
        'labeldirections' = [horizontal,vertical],
        'axesfont' = [ TIMES, ROMAN, 10 ],
        legendstyle=[location=left],
        legend=[`things vs stuff`],
        'labelfont' = [ TIMES, ROMAN, 10]):

# The following shows the labels in the exported file.
P;

# The following does not show the labels in the exported file.
plots:-display( P, 'axesfont' = [ TIMES, ROMAN, 10 ],
                'labelfont' = [ TIMES, ROMAN, 10] );

# The following shows the labels in the exported file.
plots:-display( P, 'axesfont' = [ TIMES, ROMAN, 10 ] );

Looks like Patrick got lucky, and found the buggy combination.

acer

`Maximize` is using `fdiff` to try and compute the gradient of the objective function (since no separate gradient procedure was supplied). And fdiff can raise Digits to 26 (internally, temporarily).

You can observe what happens by looking at the printout from running the code below. (Digits reverts to 15 whenever control returns from internal fdiff calls).

restart:

max2lfGn := proc (X::Vector) local a, d, e, f, H, m, s;
e := 10.^(1-floor((1/2)*Digits));
a:= Statistics:-Mean(X);
d := Statistics:-StandardDeviation(X);
f := (m, s, H) -> lfGn(X, m, s*s, H);
Optimization:-Maximize('f'(a, s, H), s = e .. 10000*d, H = e .. 1-e,
       initialpoint = {H = .75, s = d}, optimalitytolerance = 0.01, method = 'modifiednewton')
end proc:

lfGn:=proc(a,b,c,d)
  printf("%ld\n",Digits);
  -b-c+d^2;
end proc:

trace(fdiff):

max2lfGn( Vector([1,2,3,4,5,6,7]) );

Things would be faster if you ensured that the objective evaluated under evalhf, which is not the case for your `f` which relies on lexical scoping to get at lfGn and X. That might not be possible, of course, depending on what is in your lfGn procedure.

You could force Digits to be set lower *inside* the body of your lfGN procedure itself. But you may then have to worry that fdiff could sometimes fail, if it didn't attain an accuracy it expected according to the value of Digits it saw *inbound*! So you might want to craft your own gradient procedure as well, which would itself call fdiff. (Maybe see here, for some notes on that.) You first might try it without a hand crafted gradient, just be setting Digits=15 or so, right at the start of your lfGn.

acer

You had an `end` for your if-then but not for your do-loop.

Try using end if and end do instead of just end (or indent) to avoid such mix-ups.

acer