acer

33188 Reputation

29 Badges

20 years, 206 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Carl Love I did it that way so that the name UpdatingRead wouldn't itself appear in the Variable Manager, if the code was just pasted into a worksheet. It wasn't the only way to get that effect, but served that purpose.

Of course, one could instead savelib just the assigned UpdatingRead procedure to a .mla archive located in libname in future sessions.

I wanted the OP to see the addition of his .m file's missing names to the Variable Manager, using the posted Answer code in the simplest way, without the effect of the Variable Manager being cluttered by the name UpDatingRead itself.

Who knows... maybe someone will fix assign so that it too updates the Variable Manager. In which case the behaviour would be the same either way. It was intended only as the most minor aspect, and is not a big deal.

Dr. Subramanian, your line of questioning is valid, but please note that Roman's Post is about exact rational data. And your line of questioning above is about floating-point data.

You are right, Maple 2015's LinearAlgebra:-LinearSolve is missing the ability to do direct sparse (LU based) real floating-point linear solving at working precision Digits>15. That's a regression and I've submitted a bug report. It can however still do indirect, iterated sparse real floating-point linear solving using functions from the chapter F11 of the NAG library.

I mentioned above (true at the time I wrote it) that LinearSolve was using NAG routine f01brf to LU factorize a real float Matrix and routine f04axf to do the subsequent solving for a given real float RHS. In Maple 2015 this is done at hardware working precision using functions (bundled in Maple in an external shared library) from UMFPACK.

You asked above about how to re-use the real floating-point sparse LU factorization. The following sets up and uses functions from UMFPACK to accomplish that. It works in my 64bit Linux version of Maple 2015 (and seems to also work in the 64bit Linux versions Maple 18.02 back as far as, hmm, 15.01). It also seems to run ok on my 64bit Maple2015.1 on Windows 7.

I've coded it with three stages: LU factorization, linear solving for one or multiple RHSs, and freeing of the LU factorization data. The idea -- which I believe is what you are looking for -- is the ability to factorize the LHS Matrix just once and then to do several separate linear-solving steps using different RHSs. Indeed that special use pattern is the only reason to use this code. Otherwise one would simply call LinearSolve.

One must be careful with the so-called prefactored data. It's stored internally in memory, and is not explicitly assigned to any Matrix or other structure. (The module code below simply saves a memory reference. Using an invalid memory reference would likely crash.)

The code below consists of 1D Maple Notation input, interspersed with prettyprinted plaintext output.

SparseLU:=module()
   option package;
   export LU,Solve,Free;
   local extlib, Anz, NumericA, numrows;
   LU:=proc(M::Matrix(:-storage=:-sparse,:-datatype=:-float[8],
                      :-order=:-Fortran_order))
      local extfun, numcols;
      if extlib='extlib' then
         extlib:=ExternalCalling:-ExternalLibraryName("linalg",':-HWFloat');
      end if;
      extfun:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_MatFactor',
                                              extlib);
      (numrows,numcols):=op(1,M);
      (Anz,NumericA):=extfun(numrows,numcols,M);
      NULL;
   end proc:
   Solve:=proc(numrows::posint, V::{Matrix,Vector})
      local B,extfun,res;
      if type(V,'Matrix'(':-storage'=':-sparse',':-datatype'=':-float[8]',
                         ':-order'=':-Fortran_order')) then
         B:=V;
      else
         B:=Matrix(V,':-storage'=':-sparse',':-datatype'=':-float[8]',
                   ':-order'=':-Fortran_order');
      end if;
      if extlib='extlib' then
         extlib:=ExternalCalling:-ExternalLibraryName("linalg",':-HWFloat');
      end if;
      extfun:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_MatMatSolve',
                                              extlib);
      if not(Anz::posint and NumericA::posint and numrows::posint) then
          error "invalid factored data";
      end if;
      res:=extfun(numrows,op([1,2],B),Anz,NumericA,B);
      res;
   end proc:
   Free:=proc()
      local extfun;
      if not NumericA::posint then
         error "nothing valid to free";
      end if;
      if extlib='extlib' then
         extlib:=ExternalCalling:-ExternalLibraryName("linalg",':-HWFloat');
      end if;
      extfun:=ExternalCalling:-DefineExternal(':-hw_SpUMFPACK_FreeNumeric',
                                              extlib);
      extfun(NumericA);
      NumericA:=-1;
      NULL;
   end proc:
end module:
 

MM:=Matrix([[0.,0.,-25.],[-53.,-7.,0.],[0.,-70.,0.]],
           ':-storage'=':-sparse',':-datatype'=':-float[8]');
                               [ 0.      0.     -25.]
                               [                    ]
                         MM := [-53.    -7.      0. ]
                               [                    ]
                               [ 0.     -70.     0. ]

 
with(SparseLU);

                               [Free, LU, Solve]


# precompute the sparse LU factorization
LU(MM);
 
# solve MM.ans=VV with VV a Vector
VV:=Vector([1,1,1]);

                                         [1]
                                         [ ]
                                   VV := [1]
                                         [ ]
                                         [1]

ans:=Solve(3,VV);

                                [-0.0169811320754717]
                                [                   ]
                         ans := [-0.0142857142857143]
                                [                   ]
                                [-0.0400000000000000]

MM,VV; # untouched, ok

                          [ 0.      0.     -25.]  [1]
                          [                    ]  [ ]
                          [-53.    -7.      0. ], [1]
                          [                    ]  [ ]
                          [ 0.     -70.     0. ]  [1]

# verify by computing forward error
LinearAlgebra:-Norm(MM.ans - Matrix(VV));

                                      0.

# solve MM.ans=VV with VV a Matrix
VV:=LinearAlgebra:-RandomMatrix(3,2,':-datatype'=':-float[8]');

                                    [-53.    40.]
                                    [           ]
                              VV := [21.     97.]
                                    [           ]
                                    [-25.    43.]

ans:=Solve(3,VV);

                      [-0.443396226415094    -1.74905660377358 ]
                      [                                        ]
               ans := [0.357142857142857     -0.614285714285714]
                      [                                        ]
                      [ 2.12000000000000     -1.60000000000000 ]

MM,VV; # untouched, ok

                     [ 0.      0.     -25.]  [-53.    40.]
                     [                    ]  [           ]
                     [-53.    -7.      0. ], [21.     97.]
                     [                    ]  [           ]
                     [ 0.     -70.     0. ]  [-25.    43.]

# verify by computing forward error
LinearAlgebra:-Norm(MM.ans - Matrix(VV));

                                      0.

 
# free the internal LU factorization data
Free();
 
# the factorization data is gone
Solve(3,VV);
Error, (in Solve) invalid factored data

After I posted this Answer Markiyan provided a link as source of this problem:  http://math.stackexchange.com/questions/1342291/evaluate-an-integral

I had not previously seen that stackexchange post (...I only follow the "maple" tag on that site). At least one responder there used the same method as I did to arrive at Int(-1/5*sin(w)^3/w^(3/5),w=0..infinity).

I only solved that last form numerically, in my Answer above. But the cited stackexchange thread contains a hint to a way (apart from Axel's other fine suggestion to use a Laplace transform) get an exact result. The hint was that int(sin(u)*u^p,u=0..infinity) ought to be known. And Maple does know it. Now I'm sad I didn't test with the trig identity earlier, to get a difference of terms in that form. I was too happy with the form that had just the single sin call, as it would make my numeric approach with discretization easier!

Using 64bit Maple 2015.1 on Windows 7 (and it also worked when tried in my Maple 18.02),

 

restart:

P := Int(-1/5*sin(w)^3/w^(3/5),w=0..infinity);

Int(-(1/5)*sin(w)^3/w^(3/5), w = 0 .. infinity)

combine(P);

Int((1/20)*(sin(3*w)-3*sin(w))/w^(3/5), w = 0 .. infinity)

value(%);

(1/1440)*72^(4/5)*Pi^(1/2)*GAMMA(7/10)/GAMMA(4/5)-(3/40)*2^(2/5)*Pi^(1/2)*GAMMA(7/10)/GAMMA(4/5)

simplify(%);

(1/120)*Pi^(1/2)*2^(2/5)*(3^(3/5)-9)*GAMMA(7/10)/GAMMA(4/5)

 

 

Download resexact.mw

@Axel Vogt I was a little surprised that Maple didn't find an answer for int(-1/5*sin(w)^3/w^(3/5),w=0..infinity) on its own. Good job with the Laplace transform.

I got the numeric result by discretizing the oscillatory integrand and then using evalf/Sum to compute with its Levin's u-transform to accelerate convergence. Axel, you and I have discussed that approach before. In this case the period was constant (2*Pi), but for varying period NextZero may be used. But it seems this approach sometimes needs a better implementation of that acceleration algorithm than what's currently in evalf/Sum.

 

@Preben Alsholm This sounds like it fits the pattern of a known regression for 1D Maple Notation input in the Standard GUI of Maple 2015.0 which was fixed in Maple 2015.1. Namely, in a single execution group, an error gets emitted if there are multiple prompts and some aren't separated by a statement terminator. As you say, normally that is not supposed to matter.

In Maple 2015.1 it's ok. In Maple 18.02 and earlier it was ok.

@Markiyan Hirnyk It did make a difference. Instead of having to waffle around with calls to both eliminate and solve, the formula for Y as a function of only X can be done using just eliminate. I was also addressing Carl's concerns about lexicographic dependency and luck.

And now, negate theta. X stays the same and Y flips sign.

@Markiyan Hirnyk Why not eliminate both Y and theta? The result include formulas for Y in terms of X alone, and for theta in terms of X alone (and restrictions on X, which happen to be null).

restart;
eq1:=convert( X=cos(theta) + 0.8e-1*cos(3.*theta), rational):
eq2:=convert( Y=-sin(theta)+ 0.8e-1*sin(3.*theta), rational):
sols:=[eliminate({eq1,eq2},{Y,theta})]:
seq(sols[i][2], i=1..nops(sols));  # restrictions on X (NULL)
S:=seq(eval(Y,sols[i][1]),i=1..nops(sols)):
plot(S[1],X=-1-0.8e-1..1+0.8e-1,color=red);

Now, you may wish to try and show programmatically that only S[1] of the three results in sol is real. And then you might justify the negation Y=-S[1].

@Carl Love It might also be worth noting in this context, to get the x-axes alignment effect that the aligncolumns option of plots:-display provides,

_PLOTARRAY( Matrix(2,1, [P1,P2]), _ALIGNCOLUMNS(true) );

@Carl Love I checked, and see that I submitted a bug report (SCR) on the setattribute and floats problem in the Standard GUI in July 2010.

(I have an idea that it might have regressed between Maple 11 and Maple 12, but I haven't double-checked that.)

@rlewis You might want to ensure that you're using a completely new name for the .gif target file (in case the unrebooted machine is confused about the old 0 byte file -- possibl even an orphaned kernel?).

Note that for "large" 3D animations the whole export process can use a lot of memory and a lot of time. It gets more severe as the number of frames and the size of the plotted structures grow. 100 frames can be "large" in the 3D case. Also, using plottools:-sphere rather than pointplot3d symbol=solidsphere incurs more cost.

The Maple Standard GUI can also leak memory, when rendering involved 3D plots. This can hamper 3D plot animation export which itself can use a lot of resources. It can help to export from a GUI session in which you haven't otherwise display the animation inline. Or you could run the code the generates the export fully programmatically in another interface (such as the CLI Commandline Interface).

If it happens for a modest number of frames, with an entirely new filename, using the symbol=solidsphere way, then you might check that the OS is not running any errant mserver processes. Ie, if all instances of maple owned by you are shut down and some mserver process owned by you is still running, then that process could be killed (or the host restarted).

@Karen If I knew this then I'd forgotten it.

It appears to also work on ListBox entries, component tooltips, and DataTable row and column names.

And it also seems to work with colors specified in hex format.

Using 64bit Maple 2015.0 on Linux or 64bit Maple 2015.1 on Windows 7,

Download dtfun.mw

 

@Fabian The condition -u(f,g)+u(h,e)>0 is not implied by u(anything,anything) being positive.

However,

restart:

z:=((-u(f,g)+u(h,e))^(1/a))^a:

simplify(z,symbolic);

                             -u(f, g) + u(h, e)

You might also note that in some cases you can restrict the simplification (if you don't want too much of it to be done). Eg,

restart:

z:=((-u(f,g)+u(h,e))^(1/a))^a+sin(x)^2+cos(x)^2-1:

simplify(z,power,symbolic);

                                            2         2    
                 -u(f, g) + u(h, e) + sin(x)  + cos(x)  - 1


simplify(z,symbolic);

                             -u(f, g) + u(h, e)

@rlewis In my 64bit Maple 12.02 on Windows 7 it takes 1 minute for the animate call below to compute. Then, when I use commands to programmatically export to a .gif it takes 20 sec before a 0 KB file appears and then another 40 sec before the full 2.4 MB file appears.

I used plotsetup to change the plot output device to gif. That way I didn't have to use any right-click menus to export.

A few things to note about Maple 12: the plots:-display command does not applying directional lighting or surface glossiness by default (ie. opposite of Maple 2015). So I specify those as options below.

I find that the solidcircle symbol of a 3D point plot looks clunky and rough at large symbolsize. So I use plottools:-sphere below instead. Adjust its grid option to make it more or less smooth.

restart;
with(plots): with(plottools):
f1:=cos(t)/4: f2:=sin(t)/5: f3:=sin(2*t+1)/6:
pts:=[[0,0,f1],[1+f2,1,1],[-1,0,f3]]:
ptsL:=[pts[],pts[1]]:
c:=[blue,red,green]:
p1:=seq(sphere(eval(ptsL[i],t=0),0.07,color=c[i],
               grid=[30,30],style=patchnogrid),
        i=1..3):
p2:=pointplot3d(eval(ptsL,t=0),connect,thickness=3,color=black):
display(p1,p2,scaling=constrained,lightmodel=Light4,glossiness=1.0);
p:=tt->display(seq(sphere(eval(ptsL[i],t=tt), 0.07,
                          color=c[i], grid=[30,30],
                          style=patchnogrid, lightmodel=Light4),
                   i=1..3),
               pointplot3d(eval(ptsL,t=tt),
                           connect, thickness=3, color=black),
               scaling=constrained,lightmodel=Light4,glossiness=0.0):
A:=animate(p,[t],t=0..50,frames=100):
fn:=cat(kernelopts(homedir),"/rl01.gif"):
plotsetup(gif,plotoutput=fn);
A; # This should now create the .gif file, eventually.

Running the above in 64bit Maple 12.02 on Windows 7 gave me this .gif, eventually.

 

@rlewis Using Preben's code for Maple 12, and running it in Maple 12.02, I right-clicked and chose "Export" as the bottom choice in the popped-up context-menu. I selected the GIF item from the submenu of that "Export" item. I gave it a file location, etc. It takes a few moments for Maple to finish writing out the .gif file and close the write operation.

Here's what I got when it finished writing out the file:

@Carl Love I don't want to confuse Robert, but FYI it is possible in Maple 2015 to write code which embeds a PlotComponent which contains the plot-sequence-animation and which plays it automatically (ie. fully programmatically).

Robert, Carl's suggestions 1) or 2) are very much the usual way, and are what you need to do in Maple 12.

Also, Robert, right-click context-menu export of the visible plot (animation) should allow you to produce an animated .gif file, if you want.

First 344 345 346 347 348 349 350 Last Page 346 of 607