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

One reason for the omission is if the routines were added after Mark 7 (of the NAG C Library). The Maple NAG package only supports functions already present at the time of the NAG Mark 7 C Library.

Sometimes it is possible to use Maple's so-called "wrapperless" external-calling mechanism to connect with arbitrary C/Fortran functions in external shared libraries. See here for an example of that.

But sometimes a hand-crafted "wrapper" function has to be written (in C, say) so as to bridge the two. This can be necessary if the types of the external target function aren't just the simple types supported by the wrapperless external-calling mechanism. The hand-built wrapper then has to be compiled into its own shared libarary (which is linked against the NAG library).

It can sometimes also be necessary to construct a wrapper by hand if some struct or other object has to be created with one call, and then somehow passed to another call. In the case of the NAG functions, a "NAGError" or a "NAG_Comm" object might be in this class, as would be any pointer to a C struct. But maybe "NAGError" is ok, if using a default form.

These issues might not be relevant for the `comm` and `icomm` arrays formed for an initial call to f12aac and which are then passed to f12abc. And then some intermediary results from f12abc must also be passed to f12acc. Sometimes this can be mitigated by merging the calls to, say, all of f12a{a,b,c}c  external NAG functions into a single wrapper. (Sometimes that implies, eg.for F11 routines, that preconditioning must always be done for each RHS of a linear system solving computation. And that's inefficient if one wishes to repeat, at various times, with multiple RHSs. But maybe there's no such functionality penalty here, for merging access to these three routines.)

Anyway, I've been wanting to get things working for at least ARPACK (on which these and several other NAG F12 functions are based). But I haven't yet begun even to figure out what special considerations, if any, are necessary (or what non-Maple-y objects, if any, might need to get passed around.)

If you're in a big hurry, then I likely can't help much. But if you can wait a while... I might be able to help.

acer

Yes, as far as I know one has to override each dimension separately. So forcing use of new unit `DU` in place of `meter` for the dimension of length will not make DU^2 replace m^2 for the dimension of area. And so on, and so on. I don't know of a mechanism to easily replace the base unit throughout all dimensions. I guess that one could write code to generate all the new form, say for everything that Units:-GetDimensions() returns. Several "gotchas" await, for doing that. Maybe I can ask a simpler question: Do you want force to stay as Newtons (ie. the unit of N in place of kg*m/s^2) or do you want it to be scaled and represented in terms of kg*DU/TU^2?

For your information, Maple 15 has a new command, Units:-UseUnit, which can be used instead of the pair Units:-AddSystem and Units:-UseSystem in your way. It's not quite as clean and terse as it might be, since the command doesn't seem to accept multiple replacements at once. But it can be mapped.

with(Units):

AddUnit(DistanceUnit, prefix = SI, abbreviation = DU,
               context = SI, conversion = 6378.145*m);
AddUnit(TimeUnit, prefix = SI, abbreviation = TU,
               context = SI, conversion = 806.8118744*s);

#AddSystem(mysys, GetSystem(SI),
#          DU,TU,DU/TU,DU^3/TU^2,DU^2,DU^3,DU^4,TU^2,TU^3,TU^4);

#AddSystem(mysys, GetSystem(SI),
#          DU,DU^2,DU^3,TU,1/TU,DU/TU,DU/TU^2,DU^2/TU^2);

#UseSystem(mysys);

seq(UseUnit(x),
    x in [DU,DU^2,DU^3,TU,1/TU,DU/TU,DU/TU^2,DU^2/TU^2]):

with(Units:-Standard):

sqrt(2*Unit(DU^3)/Unit(TU^2)/Unit(DU));

                                            /DU\
                            1.414213562 Unit|--|
                                            \TU/

combine(Unit(N)/Unit(kg),units);

                                            /DU \
                            102.0587335 Unit|---|
                                            |  2|
                                            \TU /

-3.0*Unit(TU^2/DU^3)*Unit(DU^2/TU^2)*Unit(DU);

                                -3.000000000

It's hard to tell how you obtained that unsimplified output, with several uncombined units. As you can see above, from the last example, it does seem to simplify when entered under the Units:-Standard environment. You could also try applying combine(...,units) on it, or even just `simplify`.

This site is severely messed up and broken at present, with respect to inlining of uploaded images and worksheets. Sorry, if the inlining got worse when your post was branched.

acer

Inside the equation, on the right hand side, change p[x]+1/2...  to p[x]*Unit(bar)+1/2...

(Or some other compatible unit of pressure)

acer

If you need to do things with G(x) other than plotting, you can get the G(x) (and/or derivatives listed) by using the output=listprocedure option for dsolve.

nn:=10;
Nt:=0.5;
NB:=2.5;
Le:=2;
Pr:=2;
b:=6;

sol := dsolve({diff(G(x), x, x)+Le*f(x)*(diff(G(x), x))+(Nt/NB)*diff(T(x), x, x) = 0,
               diff(f(x), x, x, x)+f(x)*(diff(f(x), x, x))-(2*nn)/(nn+1)*(diff(f(x), x))^2 = 0,
               diff(T(x), x, x)/Pr+f(x)*(diff(T(x), x))+NB*(diff(T(x), x))*(diff(G(x), x))+
                  Nt*(diff(T(x), x))^2=0,
               G(0) = 1, G(b) = 0, T(0) = 1, T(b) = 0, f(0) = 0,
               (D(f))(0) = 1, (D(f))(b) = 0}, numeric, output=listprocedure);

useG:=eval(G(x),sol):

plot(useG, 0..b);

acer

This should optionally bypass the pop-up dialogue altogether. (So it doesn't offer the choice of a mix of some variables' ranges specified and some queried via pop-up. It's all or none.)

I used Maple 15.01 for this. The construction below will very likely not succeed in any older release that has differences in the code of the Explore:-explore internal routine. But once built and savelib'd, it would likely work "as well" in the past few releases. No promises. YMMV.

restart:
# Let's do "live surgery" on `Explore`.

kernelopts(opaquemodules=false):

inertorig:=ToInert(eval(Explore:-explore)):

inertfoo:=ToInert(
  proc(e::uneval,{variables::list(name=range(numeric)):=NoValue,
       wrapevalf::truefalse:=false})
  local v,z2,z3,z4,r,z5,z6,z7,z8,z9,z10,z11,
        z12,z13,z14,z15,z16,z17,z18,z19,uf;
  if variables=NoValue then
    BOO;
  else
    v:=map(lhs,variables);
    r:=map(rhs,variables);
    r:=[seq(op(1,x),x in r),seq(op(2,x),x in r),seq(false,x in r)];
    uf:=wrapevalf;
  end if;
end proc):

newexplore:=FromInert(subsop(1=op(1,inertfoo),
  [5,1]=NULL,[5,2]=NULL,[5,3]=NULL,
  [5,4]=op([5,1],subsop([5,1,1,2,1]=op([5,1..11],inertorig),inertfoo)),
  [5,5]=NULL,[5,6]=NULL,[5,7]=NULL,[5,8]=NULL,[5,9]=NULL,
  [5,10]=NULL,[5,11]=NULL,inertorig)):

# That's it. We could use the `savelib` command to save `newexplore`
# to a personal .mla Library archive (made with `LibraryTools:-Create`).

# This is only the most rudimentary testing.

newexplore(plot(sin(x/(a+1))+b,x=-6..6,view=[-6..6,-4..4]),
           variables=[a=1.0..4.0,b=-3.0..3.0]);

newexplore(plot(sin(x/(a+1))+b,x=-6..6,view=[-6..6,-4..4]));

acer

Sorry if I misunderstand, but how about just,

dsolve({diff(y(t), t, t)+y(t) = 0, y(0) = 0, (D(y))(0) = 1}, numeric,
       events = [[diff(y(t), t), halt]]);

Or are you worried that y(t) might level out and stay flat or just have an inflexion but not an extremum (ie. diff(y(t),t)=0 at some point, but doesn't actually change sign)?

acer

You can do the simultaneous animating using (Plot and Button) Components.

simul_animation.mw

acer

One way is to compute the LU decomposition (with output=R) and then deduce from the entries of L. See here for a routine `ElemDecomp` for doing that kind of thing.

Easier still might be to just write your own RREF routine, and inject into it calls to a procedure that forms the elementary Matrices (as it proceeds with the reduction) and stores them as a byproduct.

But if you really are just interested in the 2x2 case (it isn't clear to me, in reading your post, sorry) then you can just hard code the result. The result will contain less members depending on whether A[1,2] or A[2,1] are already zero or whether A[1,1] or (A[2,2]*A[1,1]-A[1,2]*A[2,1])/A[1,1] is 1.

A:=Matrix([[a,b],[c,d]]);

                                      [a  b]
                                 A := [    ]
                                      [c  d]

e1,e2,e3,e4 := ElemDecomp(A);

                            [1  0]          [1      0    ]  [   b]
                            [    ]  [a  0]  [            ]  [1  -]
          e1, e2, e3, e4 := [c   ], [    ], [   d a - c b], [   a]
                            [-  1]  [0  1]  [0  ---------]  [    ]
                            [a   ]          [       a    ]  [0  1]

LinearAlgebra:-Norm(e1.e2.e3.e4 - A);

                                      0
A;

                                   [a  b]
                                   [    ]
                                   [c  d]
map(normal,e1^(-1).%);

                               [a      b    ]
                               [            ]
                               [   d a - c b]
                               [0  ---------]
                               [       a    ]

map(normal,e2^(-1).%);

                               [       b    ]
                               [1      -    ]
                               [       a    ]
                               [            ]
                               [   d a - c b]
                               [0  ---------]
                               [       a    ]

map(normal,e3^(-1).%);
                                   [   b]
                                   [1  -]
                                   [   a]
                                   [    ]
                                   [0  1]

map(normal,e4^(-1).%);

                                   [1  0]
                                   [    ]
                                   [0  1]

acer

> A:=[1,3,5,2,7]:
> B:=[2,5]:

> A[B];
                             [3, 7]

Do you want list A flattened before extraction, or not?

> A:=[[1,23],[8,3],5,2,7]:
> B:=[2,5]:

> A[B];
                          [[8, 3], 7]

> ListTools:-Flatten(A)[B];
                            [23, 5]

acer