acer

32400 Reputation

29 Badges

19 years, 344 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@tomleslie Are floating-point implementations of the Heun function available in stock Matlab 2015b? Or does one have to resort to 3rd party code (such as by John Matthews)?

You appear to have omitted either a space or an explicit multiplication symbol at a few places.

For example, after the "2 u" in your very first line you have an opening round-bracket without any space or explicit multiplication symbol between them. This gets parsed as an instance of function application. Ie, as u(...) an application of the unassigned operator `u`. But you appear to have intended u*(...).

And then this issue continues throughout the document, as you subtitute in for `u`, etc.

acer

In case anyone was wondering about ways to reverse the order of integration, "by hand or using IntegrationTools".

Also, here too one must avoid use of typed literal `phi` and `theta` in the presence of escaped locals in the double integral returned by the `Flux` command.

restart;

with(VectorCalculus):

v1 := x^2 + y - sin(z):

v2 := x^2 + 1/y - 2*z:

v3 := y^2 + 3*x + z:

vv := VectorField(<v1, v2, v3>, 'cartesian'[x, y, z]):

G := Flux(vv, Sphere(<0, 0, 0>, r), inert);

Int(Int(-sin(phi)*(sin(phi)*cos(theta)^3*cos(phi)^2*r^3+sin(phi)*cos(theta)^2*sin(theta)*cos(phi)^2*r^3-cos(theta)^2*cos(phi)^3*r^3-sin(phi)*cos(theta)^3*r^3-r^3*sin(phi)*cos(theta)^2*sin(theta)+cos(theta)^2*cos(phi)*r^3+cos(theta)*sin(theta)*cos(phi)^2*r^2+cos(phi)^3*r^3-3*sin(phi)*cos(theta)*cos(phi)*r^2+2*r^2*cos(phi)*sin(phi)*sin(theta)+sin(phi)*cos(theta)*sin(r*cos(phi))*r-cos(theta)*sin(theta)*r^2-cos(phi)^2*r^2-cos(phi)*r^3-1)*r, phi = 0 .. Pi), theta = 0 .. 2*Pi)

 

Just as member vv wrote, this can be done "by hand or using IntegrationTools".

 

int(int(op([1,1],G), op(2,G)), op([1,2],G));

(4/3)*r^3*Pi+4*r*Pi

with(IntegrationTools):

H := GetIntegrand(G): # the inner integral

int(int(op(1,H), op(2,G)), op(2,H));

(4/3)*r^3*Pi+4*r*Pi

Int(Int(GetIntegrand(H), GetVariable(G)=GetRange(G)), GetVariable(H)=GetRange(H)):

value(%);

(4/3)*r^3*Pi+4*r*Pi

 

But here too care must be taken not to use literal `phi` and `theta`, unless the escaped local names in the integrand have been converted to globals. That's an inconvenient convenience of the `Flux` command.

 

Int(Int(GetIntegrand(H), theta=GetRange(G)), phi=GetRange(H)): # will go wrong

value(%); # wrong, due to local theta and phi in G

-2*sin(phi)*(sin(phi)*cos(theta)^3*cos(phi)^2*r^3+sin(phi)*cos(theta)^2*sin(theta)*cos(phi)^2*r^3-cos(theta)^2*cos(phi)^3*r^3-sin(phi)*cos(theta)^3*r^3-r^3*sin(phi)*cos(theta)^2*sin(theta)+cos(theta)^2*cos(phi)*r^3+cos(theta)*sin(theta)*cos(phi)^2*r^2+cos(phi)^3*r^3-3*sin(phi)*cos(theta)*cos(phi)*r^2+2*r^2*cos(phi)*sin(phi)*sin(theta)+sin(phi)*cos(theta)*sin(r*cos(phi))*r-cos(theta)*sin(theta)*r^2-cos(phi)^2*r^2-cos(phi)*r^3-1)*r*Pi^2

int(int(op([1,1],G), theta=0..2*Pi), phi=0..Pi); # wrong, due to local theta and phi in G

-2*sin(phi)*(sin(phi)*cos(theta)^3*cos(phi)^2*r^3+sin(phi)*cos(theta)^2*sin(theta)*cos(phi)^2*r^3-cos(theta)^2*cos(phi)^3*r^3-sin(phi)*cos(theta)^3*r^3-r^3*sin(phi)*cos(theta)^2*sin(theta)+cos(theta)^2*cos(phi)*r^3+cos(theta)*sin(theta)*cos(phi)^2*r^2+cos(phi)^3*r^3-3*sin(phi)*cos(theta)*cos(phi)*r^2+2*r^2*cos(phi)*sin(phi)*sin(theta)+sin(phi)*cos(theta)*sin(r*cos(phi))*r-cos(theta)*sin(theta)*r^2-cos(phi)^2*r^2-cos(phi)*r^3-1)*r*Pi^2

int(int(convert(op([1,1],G),`global`), theta=0..2*Pi), phi=0..Pi);

(4/3)*r^3*Pi+4*r*Pi

 


Download fluxfoo2.mw

Intsead of converting the names to `global` and typing in the literal name `phi` I could also have used IntegrationTools:-GetVariable.

with(IntegrationTools):
Change(op(1,G), cos(GetVariable(op(1,G)))=t, t);

@Axel Vogt The names `phi` and `theta` in G are locals. Presumably this is so to cover the case that the user had previously assigned values to either of them. This makes it awkward to use IntegrationTools:-Change, since its 2nd and 3rd arguments have to contain the appropriate names.


restart;

with(VectorCalculus):

v1 := x^2 + y - sin(z):

v2 := x^2 + 1/y - 2*z:

v3 := y^2 + 3*x + z:

vv := VectorField(<v1, v2, v3>, 'cartesian'[x, y, z]):

G := Flux(vv, Sphere(`<,>`(0, 0, 0), r), inert);

Int(Int(-sin(phi)*(sin(phi)*cos(theta)^3*cos(phi)^2*r^3+sin(phi)*cos(theta)^2*sin(theta)*cos(phi)^2*r^3-cos(theta)^2*cos(phi)^3*r^3-sin(phi)*cos(theta)^3*r^3-r^3*sin(phi)*cos(theta)^2*sin(theta)+cos(theta)^2*cos(phi)*r^3+cos(theta)*sin(theta)*cos(phi)^2*r^2+cos(phi)^3*r^3-3*sin(phi)*cos(theta)*cos(phi)*r^2+2*r^2*cos(phi)*sin(phi)*sin(theta)+sin(phi)*cos(theta)*sin(r*cos(phi))*r-cos(theta)*sin(theta)*r^2-cos(phi)^2*r^2-cos(phi)*r^3-1)*r, phi = 0 .. Pi), theta = 0 .. 2*Pi)

L:=[indets(G,And(name,Non(constant)))[]];

[phi, r, theta]

for nm in L do
  if addressof(nm) <> addressof(convert(nm,`global`)) then
    print(nm, addressof(nm), addressof(convert(nm,`global`)));
  else
    print(nm, "ok");
  end if;
end do;

phi, 18446884655745329310, 18446884655745292254

r, "ok"

theta, 18446884655745329278, 18446884655745328862

op(1,G);

Int(-sin(phi)*(sin(phi)*cos(theta)^3*cos(phi)^2*r^3+sin(phi)*cos(theta)^2*sin(theta)*cos(phi)^2*r^3-cos(theta)^2*cos(phi)^3*r^3-sin(phi)*cos(theta)^3*r^3-r^3*sin(phi)*cos(theta)^2*sin(theta)+cos(theta)^2*cos(phi)*r^3+cos(theta)*sin(theta)*cos(phi)^2*r^2+cos(phi)^3*r^3-3*sin(phi)*cos(theta)*cos(phi)*r^2+2*r^2*cos(phi)*sin(phi)*sin(theta)+sin(phi)*cos(theta)*sin(r*cos(phi))*r-cos(theta)*sin(theta)*r^2-cos(phi)^2*r^2-cos(phi)*r^3-1)*r, phi = 0 .. Pi)

IntegrationTools:-Change(op(1,G), cos(phi)=t, t);

Error, (in IntegrationTools:-Change) the dependent variable(s) must be of type 'unknown', as in f(x, ...) or f[i](x, ...) (no rule for the derivative of f). Received cos(phi)

IntegrationTools:-Change(op(1,convert(G,`global`)), cos(phi)=t, t);

-r*(Int((-t^2+1)^(1/2)*cos(theta)^3*t^2*r^3+(-t^2+1)^(1/2)*cos(theta)^2*sin(theta)*t^2*r^3-cos(theta)^2*t^3*r^3-(-t^2+1)^(1/2)*cos(theta)^3*r^3-r^3*(-t^2+1)^(1/2)*cos(theta)^2*sin(theta)+cos(theta)^2*t*r^3+cos(theta)*sin(theta)*t^2*r^2+t^3*r^3-3*(-t^2+1)^(1/2)*cos(theta)*t*r^2+2*r^2*t*(-t^2+1)^(1/2)*sin(theta)+(-t^2+1)^(1/2)*cos(theta)*sin(r*t)*r-cos(theta)*sin(theta)*r^2-t^2*r^2-t*r^3-1, t = -1 .. 1))

student[changevar](cos(phi)=t, op(1,G), t);

Error, (in student/changevar/SingleSolve) unable to solve cos(phi) = t for phi

student[changevar](cos(phi)=t, op(1,convert(G,`global`)), t);

Int(-((-t^2+1)^(1/2)*cos(theta)^3*t^2*r^3+(-t^2+1)^(1/2)*cos(theta)^2*sin(theta)*t^2*r^3-cos(theta)^2*t^3*r^3-(-t^2+1)^(1/2)*cos(theta)^3*r^3-r^3*(-t^2+1)^(1/2)*cos(theta)^2*sin(theta)+cos(theta)^2*t*r^3+cos(theta)*sin(theta)*t^2*r^2+t^3*r^3-3*(-t^2+1)^(1/2)*cos(theta)*t*r^2+2*r^2*t*(-t^2+1)^(1/2)*sin(theta)+(-t^2+1)^(1/2)*cos(theta)*sin(r*t)*r-cos(theta)*sin(theta)*r^2-t^2*r^2-t*r^3-1)*r, t = -1 .. 1)

 


Download fluxbah.mw

@EugeneKalentev By numeric quadrature,

evalf(Flux(vv, Sphere(<0, 0, 0>, 2.3), inert));
                          79.86766283

evalf(Flux(ee, Sphere(<0, 0, 0>, 2.3), inert));
                          79.86766283

@Carl Love I was pretty sure that's what you meant, sure.

However for this example it is faster to do evalhf(map(...)) than map[evalhf](...) on my 64bit Maple 2015 for Linux at least.

restart;
with(ImageTools):
img := Scale(Read(cat(kernelopts(datadir), "/tree.jpg")),10.0):

cimg:=CodeTools:-Usage( Complement(img) ):
memory used=243.56MiB, alloc change=241.09MiB, cpu time=190.00ms, real time=206.00ms, gc time=10.00ms

fimg:=CodeTools:-Usage( FitIntensity(img,1..0) ):
memory used=0.71GiB, alloc change=0.71GiB, cpu time=2.96s, real time=993.00ms, gc time=50.00ms

eimg:=CodeTools:-Usage( evalhf(map(x->1-x,img)) ):
memory used=243.53MiB, alloc change=243.53MiB, cpu time=2.90s, real time=2.21s, gc time=0ns

mimg:=CodeTools:-Usage( map[evalhf](x->1-x,img) ):
memory used=1.90GiB, alloc change=-243.53MiB, cpu time=7.08s, real time=6.67s, gc time=1.20s

One may also observe that neither map(x->1-x,img) nor evalf(map(x->1-x,img)) produce a result that Maple accepts as an image, since the result lacks the float[8] datatype. So I don't know what the OP intended -- there may even have been a typo.

And neither is as fast as FitItensity, which in turn is not as fast as Complement. Somewhat interestingly the former appears to be using multiple cores -- see the "cpu time" versus "real time". I ran the above on a quad-core machine. Perhaps with many more CPU cores I'd see FitIntensity perform more like Complement. Perhaps Complement could be multithreaded.

The Complement command computes faster than the other suggestions, as the color image size gets large.

Using 64bit Maple for Linux with Maple 2015,

restart;
with(ImageTools):
img := Scale(Read(cat(kernelopts(datadir), "/images/tree.jpg")),10.0):

CodeTools:-Usage( Complement(img) ):
memory used=243.56MiB, alloc change=241.09MiB, cpu time=200.00ms, real time=210.00ms, gc time=10.00ms

CodeTools:-Usage( FitIntensity(img,1..0) ):
memory used=0.71GiB, alloc change=0.71GiB, cpu time=2.99s, real time=988.00ms, gc time=50.00ms

CodeTools:-Usage( map[evalhf](x->1-x,img) ):
memory used=1.90GiB, alloc change=-243.53MiB, cpu time=8.05s, real time=6.74s, gc time=1.27s

Also, the Complement command leaves the alpha channel alone for an RGBA image. The other two suggestions made produce results that are less useful, IMHO. Those other two methods already seem more complicated to use, and fixing them up to only affect the first three channels of a color RGBA image would make then more complicated still.

restart;
with(ImageTools):

img := Scale(Read(cat(kernelopts(datadir), "/images/tree.jpg")),0.6):

L := [[img, Complement(img)],
     [FitIntensity(img,1..0), map[evalhf](x->1-x,img)]]:

Embed(L, exterior=none, interior=none);

Aimg:=ToRGBA(img);
Aimg[1..160,1..100,4]:=Array(1..160,1..100,0.7,datatype=float[8]):

AL := [[Aimg, Complement(Aimg)],
       [FitIntensity(Aimg,1..0), map[evalhf](x->1-x,Aimg)]]:

Embed(AL, exterior=none, interior=none);

ps. Thanks, Thomas, I must have pasted from the wrong sheet.

@ccAndrew How is that different from ImageTools:-Complement, where complementary colors are exchanged?

@vv Here is a hot fix to a running session, which I believe achieves your edit to the procedure (in Maple 2015.1, of course).

 

restart;

kernelopts(opaquemodules=false):

new:=FromInert(subsop([5,2]=_Inert_IF(_Inert_CONDPAIR(_Inert_INEQUAT(_Inert_LIST(_Inert_EXPSEQ(_Inert_LOCAL(15))), _Inert_LIST(_Inert_EXPSEQ())), _Inert_STATSEQ(_Inert_ASSIGN(_Inert_LOCAL(9), _Inert_FUNCTION(_Inert_ASSIGNEDNAME("nops", "PROC", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_LOCAL(15))))))),ToInert(eval(plottools:-`transform/object`)))):

unprotect(plottools);
plottools:-`transform/object`:=eval(new):
protect(plottools);

kernelopts(opaquemodules=true):

plots:-shadebetween(2, 1/sqrt(1-x^2), x=0..sqrt(3)/2);

kernelopts(version);

`Maple 2015.1, X86 64 LINUX, Jun 2 2015, Build ID 1048735`

 

 

Download shade_hotfix.mw

The `option` on a procedure can be seen with, say,

kernelopts(opaquemodules=false):
interface(verboseproc=3):

op(3,eval(plottools:-`transform/object`));

      Copyright (c) 1997 by Waterloo Maple Inc. All rights reserved.

print(plottools:-`transform/object`);

@vv That Library routine has option `Copyright (c) 1997 by Waterloo Maple Inc. All rights reserved.`.  Your (only slightly) edited version is missing the copyright.

You could alternately post a FromInert(...,ToInert(...)...) call that modifies the procedure in-session, without violating copyright.

For this example the lower curve may be added. It makes the lower edge appear smoother.

The code already included the curve along the upper boundary (which is why it seems a little darker). This now even closer to what plots:-shadebetween does internally.

(f,g) := 2, 1/sqrt(1-x^2):

plots:-display(plot(g,x=0..sqrt(3)/2),
   plottools:-transform(unapply([x,y+g],x,y))(plot(f-g,x=0..sqrt(3)/2,filled=true)));

@Carl Love Good eye! I had forgotten about that aspect of the calling sequence.

An empty set suffices to designate lack of inequality constraints. So in this example below the operator designates the equality constraint x+y=2 .

Optimization:-Minimize((x,y)->x, {}, {(x,y)->x+y-2}, -1..1, -2..2);          

                                    -15  [                     -15]
           [-0.111022302462515654 10   , [-0.111022302462516 10   ]]
                                         [                        ]
                                         [           2.           ]

An example showing this clearly on the help-page would be nice.

@Markiyan Hirnyk I followed the instruction for inequality constraints, as that is what your examples contain. 

That piece of documentation does not make much sense to me either, as far as equality constraints go.

But we can test that the form I used is being interpreted as an ineqaulity constraint (and for the examples you gave... that is relevant).

Optimization:-Minimize(x, {x+y=2}, x=-1..1, y=-2..2);          

                             -8                            -8
       [-0.103259778874000 10  , [x = -0.103259778874000 10  , y = 2.]]

Optimization:-Minimize(x, {x+y<=2}, x=-1..1, y=-2..2);

                           [-1., [x = -1., y = 0.]]

Optimization:-Minimize((x,y)->x, {(x,y)->x+y-2}, -1..1, -2..2);     

                                       [-1.]
                                 [-1., [   ]]
                                       [0. ]

The above indicates that the form of operator (for a given inequality constraint) that I used is being interpreted as the documentation describes for inquality constraints (converted to operator form). I do agree with you that the documentation that describes how equality constraints ought to be handled is quite unclear.

If I had an equality constraint in expression form then I might approach it as follows, using a pair of inequality constraints (that together imply the equality constraint).

Optimization:-Minimize((x,y)->x, {(x,y)->x+y-2,
                                  (x,y)->-(x+y-2)},
                        -1..1, -2..2);

                                    -15  [                     -15]
           [-0.111022302462515654 10   , [-0.111022302462516 10   ]]
                                         [                        ]
                                         [           2.           ]

So I will submit a but report that the documentation for expressing equality constraints in operator form is unclear and incomplete (missing).

@Markiyan Hirnyk I suspect that you have not properly realized that an inequality constraint has to be put into the form v(x1,x2,...,xn)≤0 before taking its left-hand-side as the expression to use inside the procedure.

First 320 321 322 323 324 325 326 Last Page 322 of 592