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

The problem seems not so much with `fsolve` as with what eq(0) and Y(0) might be.

If eq(0) is undefined for n=(1-k)/2 or some other value of `n` then the problem might be with your definition of `eq` or with your expectation.

Also, for your given value of `m`, `k`, and `h` it might be that Y(0) does not exist (because eq(0) may have no real root) in the range of say, n=0.47..0.51. I didn't look hard, to see whether this was a numeric roundoff issue or actually an existence issue.

acer

The right-click Export from the plot's context-menu seems to obey the dimensions of the plot as it appears in the worksheet. But you want programmatic generation of exported plots (which is entirely reasonable).

The programmatic exporting plot driver seems pretty busted in Standard. One thing I have discovered, though, is that the `height` and `width` options do work provided that you don't specify the `pt` literally when supplying the size in units of points. That is, supply it like 1014 instead of 1024pt.

This trick did not seem to succeed for the `axiswidth` and `axisheight` options, but maybe you can make do with just `width` and `height` for now.

For example, from a Worksheet in the Standard GUI of Maple 15.01, sending the plot to C:/P.ps

restart;

plotsetup(default):

P:=plot([seq(rho*sigma^2,rho=0.6..1.0,0.1)],
        sigma=0.0..0.5,
        legend=[seq('rho'=rho,rho=0.6..1.0,0.1)],
        legendstyle=[location=left]):

P;

MakeWidePlot := proc(p::evaln)
     local name, place, opts:
         name := cat("c:/temp/",convert(p,string),".ps"):
         opts := `landscape,width=1024,height=768,noborder`:
     plotsetup('ps', 'plotoutput'=cat(name), 'plotoptions'=opts):
     print( plots:-display( eval(p), 'axesfont' = [ TIMES, 10 ],
                            'labelfont' = [ TIMES, ROMAN, 10] ) ):
     plotsetup(default):
   end proc:

MakeWidePlot(P);

The image in the Postscript file (rotated right, from Portrait to Landscape) generated by the above looks like this (this is just conversion to jpeg that I did externally):

The noborder option seems to work. But the landscape option seems to have no effect (for that driver). I used 32bit Windows XP here.

A more complicated workaround could be to (in a procedure) write the plot to .m as a file, then as a system call invoke the commandline interface which would quietly read that .m and export the postscript file using its superior export driver.) Any takers?

acer

In the `LinearSolve` case, debugging can reveal that only the first two entries of the third column are being passed to the external solver. The third entry of the third column (which is the RHS of the augmented system) is passed as the value zero. And so if the external solver happens to use the second and third rows of `a` (instead of, say, the first and second) then that incorrect solution is found.

In other words, LinearSolve is computing the solution to this problem instead:

> restart:

> with(LinearAlgebra):

> a:=Matrix([[3100,6400,23610],[250,360,0]]);

                               [3100  6400  23610]
                          a := [                 ]
                               [ 250   360      0]

> evalf(LinearSolve(a));

                               [-17.56115702]
                               [            ]
                               [ 12.19524793] 

It looks like a coding error in LinearSolve, possibly made after computing the rank of `a` (which is 2, for your example).

I will submit a bug report against this.

acer

Try it instead with the option as,

mode = :-log

acer

It is not correct, to say that the `f` in the body of procedure `DO` should be subsituted by `a`. It should be substituted by the anonymous procedure (x1,x2)->a(x1,x2,o) because that is what is being passed as the first argument to `DO`, and `f` is the matching first formal parameter of `DO`.

As it turns out, the returned `f` in your result is not the same as the `f` used as the formal parameter for `DO`. It is easy to demonstrate this, by using the name `F` instead.

restart:

DO:=(F,x1,x2)->(sqrt(2)*x1+x2)*D[2](F)(x1,x2)-x2:

a:=(y1,y2,o)->1:

b:=(z1,z2,p)->simplify(DO((w1,w2)->a(w1,w2,p),z1,z2)):

b(X1,X2,OO):

                       (1/2)                             
      D[2](f)(X1, X2) 2      X1 + D[2](f)(X1, X2) X2 - X2

It can also be shown that this returned `f` is not the same as the global name f. So, this `f` in the result is some artefect from simplify, probably present because it was an unnamed procedure that got passed to `DO`. If we name it, then the 'substitution' is as expected.

restart:

DO:=(F,x1,x2)->(sqrt(2)*x1+x2)*D[2](F)(x1,x2)-x2:

a:=(y1,y2,o)->1:

thatproc:=(w1,w2)->a(w1,w2,p):

b:=(z1,z2,p)->simplify(DO(thatproc,z1,z2)):

b(X1,X2,OO);

                        (1/2)                                    
D[2](thatproc)(X1, X2) 2      X1 + D[2](thatproc)(X1, X2) X2 - X2

Now we can go back and corroborate that the mysterious `f` is just an "escaped local" from `simplify`. (Which is a bug. `simplify` is likely forgetting to eval the local name, to reproduce the anonymous proc in the returned result.)

> restart:

> DO:=(F,x1,x2)->(sqrt(2)*x1+x2)*D[2](F)(x1,x2)-x2:
> a:=(y1,y2,o)->1:
> b:=(z1,z2,p)->simplify(DO((w1,w2)->a(w1,w2,p),z1,z2)):

> hmmm:=b(X1,X2,OO);

                       (1/2)                             
      D[2](f)(X1, X2) 2      X1 + D[2](f)(X1, X2) X2 - X2

> convert(indets(hmmm),list);

                   [X1, X2, D[2](f)(X1, X2)]

> op(op(0,%[3]));

                               f

> type(%,`local`), addressof(%), addressof(f);

                   true, 162645992, 449044880

> ## It even evaluates to the original!

> eval( op(op(0,convert(indets(hmmm),list)[3])) );

                  (w1, w2) -> a(w1, w2, OO)

So, there it is. A local `f` in the result, which is a bug.

And the next disappointment might be that D[2] of that procedure does not produce zero.

acer

3) You can add a new Comment on an old Question (whether answered or not) or an Answer to such. It can help if you give more detail, or ask for clarification of details of an earlier Answer, etc.

2) Your Question which was a duplicate of this (and really, about the case of cp<3220) has been deleted (by me). Sorry, if that seems too heavy handed, but the site's had a few problematic re-posters. And all future readers benefit from keeping relevent answers grouped and organized, which doesn't happen if answers to the same question get spread across dupliacte posts unecessarily. The comment to it (by Markiyan, I think) was to the effect that no solution exists in the case at issue.

Whenever one posts any followup Comment, the old Question gets put at the top of the Active Conversations (a.k.a "recent") page. Experience shows that the members who respond on this site do see and notice such things, or often they get RSS feeds of all new posts, etc.

So my advice would be: if nobody answers, ping the original post by adding a new comment. Christopher2222 revived a dead Question just last week, after it had sat unanswered for about a month. People noticed, interest sparked, several people responded, and he seems to have solved his issue. I think it can work, just as well as can posting a duplicate, top Question.

best regards,
acer