acer

32632 Reputation

29 Badges

20 years, 45 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Axel Vogt Thanks, Axel. I believe that you once cited that document once before (somewhere), because at the time it got me interested in Horner's scheme for Matrix polynomial evaluation (at a point value for variable x), which I have code for somewhere... And that could lead to Estrin's form, and parallelization. But I digress.

I did not divide by x^3, and only subtracted by the x term when computing the initial minimax on the road to the result above. (But my motivation was also the Taylor's expansion...) A slight modification pretty easily got me what appears to be a more accurate double-precision performer as this shorter result,

proc(x::float)
 local t1;
 t1 := x^2;
 (.4290978931547805+(-0.2150922379775017e-1+(0.1194845227917730e-4
  +(0.1351777597551553e-6+0.1467942640582486e-8*t1)*t1)*t1)*t1)*t1*x
  /(1.287293679464341+(-.5794451431789617+0.2339492595289925e-1*t1)*t1)
  +x;
end proc

Getting some reasonably simple result that could attain, oh say, 0.6 ulps maximal error on [0..Pi/4] would be nice to show as Maple code.

acer

@Axel Vogt Thanks, Axel. I believe that you once cited that document once before (somewhere), because at the time it got me interested in Horner's scheme for Matrix polynomial evaluation (at a point value for variable x), which I have code for somewhere... And that could lead to Estrin's form, and parallelization. But I digress.

I did not divide by x^3, and only subtracted by the x term when computing the initial minimax on the road to the result above. (But my motivation was also the Taylor's expansion...) A slight modification pretty easily got me what appears to be a more accurate double-precision performer as this shorter result,

proc(x::float)
 local t1;
 t1 := x^2;
 (.4290978931547805+(-0.2150922379775017e-1+(0.1194845227917730e-4
  +(0.1351777597551553e-6+0.1467942640582486e-8*t1)*t1)*t1)*t1)*t1*x
  /(1.287293679464341+(-.5794451431789617+0.2339492595289925e-1*t1)*t1)
  +x;
end proc

Getting some reasonably simple result that could attain, oh say, 0.6 ulps maximal error on [0..Pi/4] would be nice to show as Maple code.

acer

I'm not sure whether the grid/cell structure in use is related to the differing behaviours. The attached worksheet gives another variable view, and overlays with the `outlines` option. It seems to me that for some values of the parameters there may be visual indications of a tie-in between the cell structure and the behaviours. But as Alejandro mentioned, the source code in question may be compiled and unavailable to view, so that can only be conjecture.

interpexplore.mw

acer

Perhaps this simpler example, animated, reveals what is going on.

As the two parallel lines (seen when `a` is farther away from the value 0) get closer together there is a point when something curious happens.

The distance between the computed cells in each of the two solution lines becomes greater than (something like) the distance between those two lines. And for `a` less than that critical value, the routine appears to be identifying the solution as a set of islands instead of as two lines. Points from both lines suddenly start getting identified as being within the same cell. (I'm not sure whether "cell" is the right term here. It is used in the ?implicitplot documentation, and possibly this mirrors what happens in the Arrays in the CURVES section of the generated plot structure.)

This effect can be delayed a bit, by using a huge `numpoints` and large `gridrefine` setting. But eventually, as the constant term gets smaller the curiosity occurs.

And then, as `a` continues to shrink toward 0 value, those solution islands shrink to points (and might even vanish).

plots:-animate(plots:-implicitplot,[(x-y)^2+a, x = 0 .. 1,
                                     y = 0 .. 1],
               a=-0.01..0.0, frames=100);

(This is essentially the same as Preben's animation, simplified and with a different slope.)

acer

@ThU Well, if you only ever want the shared colour coding scheme to be ZHUE (hue determined by z-value) then there is a much easier way. The stuff I wrote above is something I cooked up with more general coloring schemes in mind.

 

restart:
f1 := (x, y) -> sin((1/2)*y+(1/2)*x):

f2 := f1+1:

plots:-display(plot3d(f1, 0..10, 0..20),
               plot3d(f2, 0..10, 0..20),shading=zhue,
               orientation=[30,70,10]);

 

Download sharedzhue.mw

@ThU Well, if you only ever want the shared colour coding scheme to be ZHUE (hue determined by z-value) then there is a much easier way. The stuff I wrote above is something I cooked up with more general coloring schemes in mind.

 

restart:
f1 := (x, y) -> sin((1/2)*y+(1/2)*x):

f2 := f1+1:

plots:-display(plot3d(f1, 0..10, 0..20),
               plot3d(f2, 0..10, 0..20),shading=zhue,
               orientation=[30,70,10]);

 

Download sharedzhue.mw

@Markiyan Hirnyk Yes, I thought that was obvious.

I changed f2 from f1+2 to be f1+1, so that there was overlap. I think that the overlap illustrates better that the colors match for any given shared z-value, since with overlap more z-values (than just the singleton z=1) are shared.

And, of course, maxz and minz must get changed accordingly, if they are to be force fed for an example. I gave the appropriate maxz and minz for my example, with the different f2=f1+1.

This also relates to having a more general scheme (mine given at top being non-bug-free, I'm sure) which can extract the maximal and minimal z-values automatically.

@Markiyan Hirnyk Yes, I thought that was obvious.

I changed f2 from f1+2 to be f1+1, so that there was overlap. I think that the overlap illustrates better that the colors match for any given shared z-value, since with overlap more z-values (than just the singleton z=1) are shared.

And, of course, maxz and minz must get changed accordingly, if they are to be force fed for an example. I gave the appropriate maxz and minz for my example, with the different f2=f1+1.

This also relates to having a more general scheme (mine given at top being non-bug-free, I'm sure) which can extract the maximal and minimal z-values automatically.

@ThU I tried to show a more general mechanism, where one does not know in advance what the max and min values -- global over both f1 and f2 combined -- will be. Hence I dug it out of the initial plots. If one knows those in advance then the coding can be simpler.

And it can also be done without all the `op` and data structure manipulation. And it can look a little simpler with expressions instead of operators (...I think). Note that I didn't take any special care for efficiency before, either.

 

restart:

f1 := sin((1/2)*y+(1/2)*x):
f2 := f1 + 1:

maxz, minz := 2, -1:

p1 := plot3d(f1, x=0..10, y=0..20,
             color=[(maxz-f1(x,y))/(maxz-minz), 0.75, 0.75,
                    colortype=HSV]):
p2 := plot3d(f2, x=0..10, y=0..20,
             color=[(maxz-f2(x,y))/(maxz-minz), 0.75, 0.75,
                    colortype=HSV]):

plots:-display(p1,p2,axes=box,orientation=[25,70,10]);

restart:

f1 := sin((1/2)*y+(1/2)*x):
f2 := f1 + 1:

maxz, minz := 2, -1:

F := fun -> plot3d(fun, x=0..10, y=0..20,
                 color=[(maxz-fun)/(maxz-minz), 0.75, 0.75,
                         colortype=HSV]):

plots:-display(F(f1),F(f2),axes=box,orientation=[25,70,10]);

 

 

Download colorfun2b.mw

@ThU I tried to show a more general mechanism, where one does not know in advance what the max and min values -- global over both f1 and f2 combined -- will be. Hence I dug it out of the initial plots. If one knows those in advance then the coding can be simpler.

And it can also be done without all the `op` and data structure manipulation. And it can look a little simpler with expressions instead of operators (...I think). Note that I didn't take any special care for efficiency before, either.

 

restart:

f1 := sin((1/2)*y+(1/2)*x):
f2 := f1 + 1:

maxz, minz := 2, -1:

p1 := plot3d(f1, x=0..10, y=0..20,
             color=[(maxz-f1(x,y))/(maxz-minz), 0.75, 0.75,
                    colortype=HSV]):
p2 := plot3d(f2, x=0..10, y=0..20,
             color=[(maxz-f2(x,y))/(maxz-minz), 0.75, 0.75,
                    colortype=HSV]):

plots:-display(p1,p2,axes=box,orientation=[25,70,10]);

restart:

f1 := sin((1/2)*y+(1/2)*x):
f2 := f1 + 1:

maxz, minz := 2, -1:

F := fun -> plot3d(fun, x=0..10, y=0..20,
                 color=[(maxz-fun)/(maxz-minz), 0.75, 0.75,
                         colortype=HSV]):

plots:-display(F(f1),F(f2),axes=box,orientation=[25,70,10]);

 

 

Download colorfun2b.mw

I would say it even more strongly: it is a much worse idea to try and integrate and differentiate from an interpolated set of discrete data values that it is to instead work directly with the differential system.

So, if one has an ode system in Maple then do not first extract discrete values of the solution and then interpolate, and then integrate/differentiate/root-find/etc. It is generally a terrible idea to do so, from a numerical standpoint.

For integraton and differentiation it is likely a better idea in general to augment the ode system, and one should be able to find several posts about that on this site. For multiple root-finding of some expression in terms of the solution it is generally better to make the extra effort and use dsolve,events.

Suppose that you want to integrate the solution to some ode. Three possible ways are: 1) generate discrete data, interpolate, and integrate that interpolating result, 2) obtain a procedure from dsolve/numeric and then use evalf/Int on that, and 3) augment the ode with a new variable and a new equation (which sets that new variable equal to the derivative of what you want integrated) and call dsolve/numeric on the augmented ode.

Of these ways, 1) is by far the worst, numerically, in general.

I believe that 3) is generally better than 2), numerically, because the dsolve/numeric engine knows best how to do its own error analysis and step-size control dynamically to attain any requested tolerance. Whereas passing a piecewise interpolating procedure (as returned by dsolve/numeric) on to some numeric integrator introduces an unnecessary break in the informational flow. The evalf/Int integrator will be generally restricted in its accuracy by the fixed quality of the returned dsolve/numeric procedure.

Differentiation near the end-points of the original region of interest will be especially difficult to do effectively using an interpolating scheme from a precomputed set of data from the numerically solved ode.

acer

I would say it even more strongly: it is a much worse idea to try and integrate and differentiate from an interpolated set of discrete data values that it is to instead work directly with the differential system.

So, if one has an ode system in Maple then do not first extract discrete values of the solution and then interpolate, and then integrate/differentiate/root-find/etc. It is generally a terrible idea to do so, from a numerical standpoint.

For integraton and differentiation it is likely a better idea in general to augment the ode system, and one should be able to find several posts about that on this site. For multiple root-finding of some expression in terms of the solution it is generally better to make the extra effort and use dsolve,events.

Suppose that you want to integrate the solution to some ode. Three possible ways are: 1) generate discrete data, interpolate, and integrate that interpolating result, 2) obtain a procedure from dsolve/numeric and then use evalf/Int on that, and 3) augment the ode with a new variable and a new equation (which sets that new variable equal to the derivative of what you want integrated) and call dsolve/numeric on the augmented ode.

Of these ways, 1) is by far the worst, numerically, in general.

I believe that 3) is generally better than 2), numerically, because the dsolve/numeric engine knows best how to do its own error analysis and step-size control dynamically to attain any requested tolerance. Whereas passing a piecewise interpolating procedure (as returned by dsolve/numeric) on to some numeric integrator introduces an unnecessary break in the informational flow. The evalf/Int integrator will be generally restricted in its accuracy by the fixed quality of the returned dsolve/numeric procedure.

Differentiation near the end-points of the original region of interest will be especially difficult to do effectively using an interpolating scheme from a precomputed set of data from the numerically solved ode.

acer

@maple_matt_2 The subexpression 1/8*M[b]*(1+1/``(tan(theta))^2) does not appear explicitly in EQN. You would have had to copy the 1/8*M[b] term and the bracketed ``(1+1/``(tan(theta))^2) term seperately (or done something similar) as the blue output of EQN has a V[c] term lying between those in the product subexpression.

But, sure, entering an explicit multiplication symbol may be better practice than trying to remember whether you had instead added a space between paste-ins.

This problem seems to come up quite often. Some people configure to use 1D Maple notation for input, because of the risk here.

@maple_matt_2 The subexpression 1/8*M[b]*(1+1/``(tan(theta))^2) does not appear explicitly in EQN. You would have had to copy the 1/8*M[b] term and the bracketed ``(1+1/``(tan(theta))^2) term seperately (or done something similar) as the blue output of EQN has a V[c] term lying between those in the product subexpression.

But, sure, entering an explicit multiplication symbol may be better practice than trying to remember whether you had instead added a space between paste-ins.

This problem seems to come up quite often. Some people configure to use 1D Maple notation for input, because of the risk here.

@Christopher2222 Rebuilding the data file is not difficult.

I ran these commands, at the end of the Document (and then deleted those commands). I changed the data file name from "test.dat" to the more Windows-friendly "testdata.txt".

stuff:=DocumentTools:-GetProperty('plotter_',value):
writedata(cat(kernelopts(homedir),"/testdata.txt"),
          convert(op([1,3],[plottools:-getdata(stuff)]),listlist),
          float);

Since the above creates the data file "testdata.txt" in one's home directory/folder, I figure that would also be an OK place to put the "filenames.txt" file which is the applications master record of any data files it is supposed to handle. Therefore I added a line to the Document's hidden (code) document block to set the working directory to be cat(kernelopts(homedir),"/test.dat") so that it could find these two text files.

Here are the edited files. I removed the plot and output from the Document, so that one can see it changed as it works for the first time. The first two should be placed in whatever Maple reports as the output from running the command kernelopts(homedir).

filenames.txt

testdata.txt

DataAnalysis_edit.mw

First 402 403 404 405 406 407 408 Last Page 404 of 597