acer

32490 Reputation

29 Badges

20 years, 8 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by 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

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

B:=Matrix(2, [[e,f],[g,h]]):

A . B;

                           [a e + b g  a f + b h]
                           [                    ]
                           [c e + d g  c f + d h]

LinearAlgebra:-Multiply(A,B);

                           [a e + b g  a f + b h]
                           [                    ]
                           [c e + d g  c f + d h]

LinearAlgebra:-MatrixMatrixMultiply(A,B);

                           [a e + b g  a f + b h]
                           [                    ]
                           [c e + d g  c f + d h]

with(LinearAlgebra):

Multiply(A,B);

                           [a e + b g  a f + b h]
                           [                    ]
                           [c e + d g  c f + d h]

MatrixMatrixMultiply(A,B);

                           [a e + b g  a f + b h]
                           [                    ]
                           [c e + d g  c f + d h]

You can make any of those into an assignment to the name C, using := in the usual way.

acer

What would a*a produce?

How about somehow registering a name. You could give it an attribute, for example, as the mechanism. And this would allow you to specify the order, too. (I haven't dealt with the mixed case, but you could alter the code to do so.)

It's unclear what you'd want done with c*d, for protected c and d, or for any other such special pair. What decides which product produces the 1, and which the -1? Alphabetic order, say, seems an strange control, but what do I know? This too could be accomplished: use just one kind of attribute to mark names as special, and then test their alphabetic precedence in the code's conditionals.

Note that :-`*` is the global (usual) star operator. So that's what this version falls back to, in cases which it doesn't recognize.

restart:

m:=module() option package;
   export `*`;
   `*`:= proc()
      if nargs=2 and type(args[1],name)
        and type(args[2],name) then
         if attributes(args[1])=lefty
           and attributes(args[2])=righty then
            return 1;
         elif attributes(args[2])=lefty
           and attributes(args[1])=righty then
            return -1;
         elif attributes(args[1])=lefty
           or attributes(args[2])=lefty
           or attributes(args[1])=righty
           or attributes(args[2])=righty then
            return 'procname'(args); # or :-`*`(args), or whatever you decide
         end if;
      end if;
      return :-`*`(args);
   end proc;
end module:

with(m);
                              [*]

17*1.2;
                              20.4

a*b; # nothing special about these yet
                              a b

b*a;
                              a b

setattribute(a,lefty):
setattribute(b,righty):

a*b;
                               1

b*a;
                               -1

a*b*a; # it inherits the "usual" associativity

                               a

b*a*b;
                               -b

An interesting situation might be where you have a*a returned untouched, and you want to be able to right-multiply this later on by b. (This doesn't happen, if you decide to code it to return say 1 for input a*a.) In such cases, you might need to start ripping apart unevaluated `*` calls when they are the inputs, to figure out what it should do. Eg,

a*a;
                            *(a, a)

%*b; # should it handle this more cleverly?
                           *(a, a) b

Have you considered examining the methodology used in the source of the quaternions package on the application center?

acer

The combination of two files, xxxx.lib and xxxx.ind, for a Library archive is out of date and should be deprecated. It was an outright bug in recent vesions up to Maple 14.01 that help-pages such as ?LibraryTools,Create still described using xxx.lib in Examples. This was fixed in Maple 15 and the Online Help, I think. The recommended way (for some years and major releases now) is to create a single xxxx.mla file for the Library archive.

A file xxxx.hdb is a help database file. See ?INTERFACE_HELP or ?makehelp.

It's really not a good idea to stick new files into your Maple installation's "lib" folder. You or someone else will end up forgetting all that you've added, some day. Better is to add it to another folder altogether, and then augment `libname` in worksheets intended to use it (eg. in Startup code, if you want that hidden from view).

You can also see chapter 10 (and especially section 10.3) of the new Programming Manual. That gets into some detail, as far as the Library side aspects of code autoring goes. But it doesn't show or link so nicely to material on the related aspect of Help Page authorship.

A good start-to-finish example worksheet, clearly covering all these topics in one place with a comprehensive but more straightforward example, on these topics would be a nice addition.

I believe that the InstallerBuilder is an underappreciated part of these tasks, as it allows one to build a self-unpacking .mla which can easily be configured to unpack all the .mla, .hdb, supporting example, source, and other files right into a neat "toolbox" folder. This has the benefit that libname automatically gets the child "lib" folder from such unpacked installations. But this too need an example, as its own help-pages are a bit obscure.

acer

Maple already has CLAPACK (v.3.0 I think) bundled.

Here's an example of setting up such a function for use as a Maple procedure.

Let us know, if you encounter troubles.

If you have a supported NAG Library installed and working, then the NAG package may allow such access already (functions from its chapter F06, F07. and F08).

Who thinks that a complete set of Maple entry points for all BLAS and CLAPACK would be worth the effort? (This should be easy to develop, as the appropriate wrappers from that NAG package could probably be used as templates.)

acer

Did you try loading Units:-Standard first? Doing so should make the multiplication of units combine automatically.

It is not true that ASCII m is an SI prefix equivalent to micro. It's the prefix for milli. See the Units,Prefixes help page. That page has a note,

  ...The correct symbol for the prefix micro is the greek letter m. Because the SI does not
 give an acceptable alternative in an ASCII environment, three prefix symbols have gained
 acceptance in various fields: u, mu, and mc. Any of these prefix symbols is valid in the
 Units package.

It looks like there is a problem with using the prefix `mu`, however, which might be related to why mu doesn't work either. This might be partially corrected by introducing a new unit muF. (Replacing the `farad` unit, to include the alternate prefix, might also be possible but I'm not sure...) This doesn't help with regard to applying the prefix to other units.

 

Error, (in Units:-Unit) `muF` is not a recognized unit

Error, (in Units:-Unit) ``μF`` is not a recognized unit

Units:-AddUnit(muF,context=SI,conversion=microfarad,symbol=`μF`,
               spellings={`&muF`});

 

Download microfarad.mw

acer

> counter:=proc() global c; c:=c+1; NULL; end proc:

> c:=0:

> [seq([seq(op([i,counter()]), j = i .. 5)], i = 1 .. 5)];

    [[1, 1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3], [4, 4], [5]]

> c;
                               15

Maybe simpler,

> restart:

> counter:=proc(x) global c; c:=c+1; x; end proc:

> c:=0:

> [seq([seq(counter(i), j = i .. 5)], i = 1 .. 5)];

    [[1, 1, 1, 1, 1], [2, 2, 2, 2], [3, 3, 3], [4, 4], [5]]

> c;
                               15

Here, the argument to `counter` is returned unchanged. So whetever lies at the heart of your innermost `seq` call can just be wrapped in a call to `counter`.

acer

First 277 278 279 280 281 282 283 Last Page 279 of 337