You attached a data file, but you stil have not attached any worksheet showing your actual code.

In a few of your responses you included something like code, as text. But it was all incomplete, and with different sizes than the attached data, and very incomplete since none of it worked alone. And there were syntax errors in extracting particular columns.

Anyway, the following works for me. I used LinearAlgebra:-Dimensions to pick off the number of rows and columns of the imported data. I used square brackets to index into the Matrix -- over a range -- to extract a single column as Vector.

signalplot_ac.mw

Yes, there are problems here. For example, using Maple 2019.2 or Maple 2020.0,

restart;
assume(0 < a, a < R);

K := 2*Pi*Int(sin(phi)
              *Int(rho^2,
                   rho = sin(phi)*R - sqrt(sin(phi)^2*R^2 - R^2 + a^2)
                         ..
                         sin(phi)*R + sqrt(sin(phi)^2*R^2 - R^2 + a^2)),
              phi = arccos(a/R) .. Pi - arccos(a/R)):

value(K);

                -Pi^2*R*a^2

#
# The following call to `simplify` merely changes the radical
# in the inner limits of integration to,
#
#    ( -R^2*cos(phi)^2 + a^2 )^(1/2)
#
simplify(K, trig):
value(%);

                2*Pi^2*R*a^2

#
# The following call to `combine` merely changes the radical
# in the inner limits of integration to,
#
#    ( -2*R^2 + 4*a^2 - 2*R^2*cos(2*phi) )^(1/2)/2
#
combine(K):
value(%);

                2*Pi^2*R*a^2

Below we can seen that the inner integral is not the problem. The problem comes from the elliptictrig method which is one of those that int tries for the outer integral, and whose result -- when it does not fail -- gets preferred and returned.

One possible way to avoid the problem is to force int to use a method other than the problematic elliptictrig. This means not using the nicely typeset 2D Input integral, so I'll use plaintext 1D Maple Notation instead. Below I'll continue in the same session, with the same assumptions on a and R.

PP := int(rho^2,
          rho = sin(phi)*R - sqrt(-R^2*cos(phi)^2 + a^2)
                ..
                sin(phi)*R + sqrt(-R^2*cos(phi)^2 + a^2));

     1/3*(sin(phi)*R+(-R^2*cos(phi)^2+a^2)^(1/2))^3
    -1/3*(sin(phi)*R-(-R^2*cos(phi)^2+a^2)^(1/2))^3

P  :=  int(rho^2,
           rho = sin(phi)*R - sqrt(sin(phi)^2*R^2 - R^2 + a^2)
                 ..
                 sin(phi)*R + sqrt(sin(phi)^2*R^2 - R^2 + a^2));

     1/3*(sin(phi)*R+(sin(phi)^2*R^2-R^2+a^2)^(1/2))^3
    -1/3*(sin(phi)*R-(sin(phi)^2*R^2-R^2+a^2)^(1/2))^3

# Those are the same. The problem is with the outer integral.
simplify(PP-P);

                     0

int(P*sin(phi),phi=arccos(a/R)..Pi - arccos(a/R), method=elliptictrig);

                  1       2
                - - R Pi a 
                  2        

# The following returns unevaluated.
int(PP*sin(phi),phi=arccos(a/R)..Pi - arccos(a/R), method=elliptictrig):
op(0,%);

                    int

int(P*sin(phi),phi=arccos(a/R)..Pi - arccos(a/R), method=ftoc);

                        2
                  R Pi a 

int(PP*sin(phi),phi=arccos(a/R)..Pi - arccos(a/R), method=ftoc);

                        2
                  R Pi a 

When no method is specified the integration of ("simplified" form) P produces the wrong result from the ellitictric method, which gets selected and returned automatically. When no method is specified the integration of the (unsimplified form) PP produces the acceptable result from the ftoc method, since the elliptictrig method fails to produce anything.

Thus a workaround is to force the ftoc method for such examples that otherwise go wrong with the elliptictrig method. That may not always convenient to set up, however.

I have previously seen several similar problematic examples with that method. I'll submit a bug report.

int_oddity.mw

Here is the original nested integral, using 2D Input and with the inner integral typeset and the outer integral forcing method=ftoc. The limits of integration are as in the original, and the desired result attains.

int_ftoc.mw

My personal preference is for a cutout view, so that any relationship between the particular edges of the washers with the surfaces can be seen.

In order to keep the GUI responsive, it helps considerably to keep the frame data smaller. By using a grid=[2,49] option for the washer construction they can retain the same look as default, while taking less memory resources. The benefit is great enough that a modest number of washers can be shown together in the frames.

restart;
F:=proc(yy, ea) local hh,opts,phi,r,y; uses plots;
  opts := style=surface, color=gold, grid=[2,49];
  display(seq([
    display(plot3d([[r*cos(phi),r*sin(phi),y],
                    [r*cos(phi),r*sin(phi),y+h]],
                   r=y^2..sqrt(8*y), phi=0..ea, opts),
            plot3d([[y^2*cos(phi),y^2*sin(phi),hh],
                   [sqrt(8*y)*cos(phi),sqrt(8*y)*sin(phi),hh]],
                   hh=y..y+h, phi=0..ea, opts)),
    display(polygonplot3d([[y^2,0,y],[y^2,0,y+h],
                           [sqrt(8*y),0,y+h],[sqrt(8*y),0,y]]),
            polygonplot3d([[y^2*cos(ea),y^2*sin(ea),y],
                           [y^2*cos(ea),y^2*sin(ea),y+h],
                           [sqrt(8*y)*cos(ea),sqrt(8*y)*sin(ea),y+h],
                           [sqrt(8*y)*cos(ea),sqrt(8*y)*sin(ea),y]]),
      transparency=0.5)][], y=0..yy, h));
end proc:
numframes:=15: h:=2.0/(numframes): endangle:=4*Pi/3:
P:=plot3d([[r*cos(phi),r*sin(phi),r^2/8],[r*cos(phi),r*sin(phi),r^(1/2)]],
          r=0..4, phi=0..endangle,
          style=surface, color=["Green","Blue"], transparency=0.75):
plots:-animate(F,[y, endangle], y=0..2-h, frames=numframes,
               background=P, paraminfo=false,
               lightmodel=light2, glossiness=1.0, scaling=constrained,
               axes=normal, labels=[z,x,y], orientation=[-60,70]);

washers_ac.mw

Naturally, it can be done similarly for shells.

@mikerostructure Yes, I have Maple 2020.0, and I use right-click to export 3D plots on a regular basis.

I usually export to .png format. If your problem is with .eps export specifically then this would a good time to mention it. (For me, .eps export of your simple example does make my machine's CPU ramp up for a couple of seconds -- the export file is not small. But it doesn't crash.)

Another possible point of interest may be whether you have a particular language setting in Maple's GUI.

In your followup comments you've started mentioning a problematic file. Is there a reason that you cannot upload and attach that particular file to a Comment here, using the green up-arrow of the Mapleprimes editor? If the problem occurs when you re-execute it from scratch then perhaps you could even "Edit/Remove All Output" within it, for an easier, smaller upload.

This is numerical error, yes. The particular behavior depends on the working precision.

Your Document has a line that assigns to the name Digit . If you changes that to be Digits instead then you can easily adjust the working precision and see what happens when it is run with Digits=15,16,17,18, etc.

If this expression is assigned to name expr, say, then you can compare with expand(expr) and combine(expand(expr)). Compare with the values that occur in the numerator of your original expression.

Note that by default this plotting call will try and utilize evalhf to do the plotting in hardware double precision. That switches to so-called software floating-point when Digits>15.

This particular plotting command will also switch from evalhf  to sofware floats as a fallback if it detects a Float(undefined) result. There is a small range from about t=14880.0 to t=14902.7 in which the expression evaluates to 0.0 because roundoff error has occured but yet not underflow/overflow. For t above that range it will fall back to software float mode (and the curve again appears even at 0.45). This attachment demonstrates this behaviour 1_ac2.mw

The particular fashion in which numerical error occurs for this example is different for forced software floats at Digits<15 versus hardware double precision. There are ways to force software floats for Digits<=15 . (see attachment)

(Inlining documents to Mapleprimes seems to be broken right now.) Here is an attachment that demonstrates that forcing of software floats at lower Digits : Download 1_ac.mw

If your plot was constructed using the plot command then all you have to do is look at the product documentation, in order to see how to specify a particular range for the independent variable.

For example the Help page for the plot command shows several examples of this basic functionality of specifying the range over which to plot.

Or, was your plot generated in some other way than typing in a command? Did you use an item from the context-menu?

Just for fun, making it look closed.

I did this using Maple 16. If your version is older and this doesn't work then tell us what version you use.

F:=plottools[transform]((x,y)->[x,0,y+x^2/8])(plot(r^(1/2)-r^2/8,r=0..4,
                                                   filled,transparency=0.0)):
P := A -> plots:-display(orientation=[-45,70],
            plot3d([[r,theta,r^(1/2)],[r,theta,r^2/8]],
                   r=0..4, theta=0..A, coords=cylindrical,
                   orientation=[-50,80], lightmodel=light1,
                   grid=[40,max(2,ceil(evalf(49*A/(2*Pi))))],
                   scaling=constrained, axes=normal, style=patch),
            `if`(A<6.28,[F,plottools[rotate](F,A,[[0,0,0],[0,0,1]])][],NULL)):
plots[animate](P, [A], A=0..2*Pi, frames=17, paraminfo=false);

Download solid_rev.mw

@Mike Graber You wrote this:

        top level.

        A:=7:

        B:=A:

        A:=9:

        C:=A:

        At this point in the program, both B and C have the value 9. 

That code, by itself, will result in B having 7 as its assigned value (ie. B will evaluate to 7), and not 9.

Perhaps you've described what you want to see happen, and not what you've seen happen. If that is what you want then I think you could have expressed it more clearly.

In that case what you want is for B to get assigned a reference to A by name (and not its prior value), so that it subsequently evaluates to whatever value A subsequently evaluates (in an ongoing manner). Here is an example of that,

A := 7;

           A := 7

B := 'A';

           B := A

A := 9;

           A := 9

B;

              9

Your use of the phrase "by address" could be better replaced with "by reference".

To access or refer to a variable by its name -- rather than by its value -- use single right-quotes. Those are sometimes called unevaluation quotes in Maple context, as they delay the evaluation of the name.

See Chapter 3 of the Programming Manual (esp. Sections 3.3 and 3.4) as one place to get some background.

See also the Help pages for Topics procedure and evaluation for some additional information on specifc behavior.

These are not dumb questions. But you really ought to provide the full details of the preliminary computations.

Use the green up-arrow in the Mapleprimes editor to upload and attach your actual worksheet.

Use colon-equals := to make an assignment to a name. Your original image (now edited in the original Question) showed phi = 10 rather than an assignment statement.

By using only equals = you create an equation, which does not assign any value to the name.

Ie ,

phi:=10

versus

phi=10

The second of those -- by itself--  has not affect on later use of phi.

You can, however, utilize the equation phi=10 within an `eval` call.

As for `assume` it is tricky. Best if you show your problem case in full. My guess is that you are running into the fact that assumptions do not prevent any assignment, and that placing new assumptions (using `assume`) and assignment will wipe previous `assume` effects. See the `additionally` command. Or utilize `assuming` instead.

You already knew about assignment statements, and I suspect that you also knew about ways to concatenate names. You've put these actions together.

Note that it is only practical in a modest number, and since the generated names are always global it's far less useful as a programming technique.

Also, each global new name will get checked against the assigned names stored in all .mla and .m files in libname (once at least, up to first evaluation). That seems to scale roughly linearly, IIRC. I once saw an example with a 15sec wait in the computation where the kernel was doing nothing but i/o overhead.

So it's sorta, kinda ok as an off-the-cuff thing at the top-level, but in my opinion doesn't belong in a group of solid programming techniques for serious production line work.

You might consider this (although nothing is fool-proof),

  select(u->is(u,real), evalc~(allvalues~(coef7))) assuming alpha[1,2]::real, alpha[2,6]::real;

or, since it is only the two assumed names as dependent names,

  select(u->is(u,real), evalc~(allvalues~(coef7))) assuming real;

I declared z as local outside the do-loop as, alas, the 2D Input parser doesn't accept your original.

Download G2_ac.mw

You might consider editing your procedures in 1D plaintext Maple Notation (in Execution Groups, or Code Edit Regions).

ps. I am curious, is your task going to go anywhere near the LinearAlgebra[Generic] subpackage? (I know nothing of "Domineering".)

You haven't stated why LinearAlgebra:-LinearSolve is not adequate for your computations. You haven't mentioned whether you are doing purely floating-point computations, and if so whether those are purely real.

I would not be surprised if you needed some large examples (or a great many of them) before the LinearSolve overhead became a bigger bottleneck than anything else in your code -- as a tentative judgement based from general experience with member's code.

But, below is an example that uses the floating-point real tridiagonal linear solver from LAPACK.

This uses two functions, one to do the LU decomposition and another to do the back/forward substitutions for one or more RHSs.

Hopefully I got the integer array's widths correct -- I might have forgetten some detail of the linkage model in use.

As I mention below, this has a little extra overhead (probably avoidable, using other access means).

This was run on Linux. To run on MS-Windows you'd have to changed the name of the external shared library, from "libmkl_rt.so" to the appropriate .dll. And even then it might run amok of MS-Windows awkwardness.

Download DGTTRF.mw

There are a few things that can be said.

1) No, it is not necessary to use active int for this. In fact doing so can be unnecessarily very expensive, since it can cause the integral to be attempted symbolically for each new call to the objective procedure (where each might fail, after a while). See first example in the attached.

There can sometimes be a mixup for Mimimize when called with procedure (rather than expression) form for the objective. It happens more rarely then it did long ago. It's caused mostly by some muddle about generating a procedure for evaluating the derivative numerically, via automatic differentiation. My first example does a few things to avoid that.

2) Your second bullet point sounds like it might be an evaluation mix-up. See the last example in the attachment with proc4, which may be what you were looking for here. The point is to keep straight that the K,r,x in globally defined df are themselves global (whereas the K,r in proc4 are procedure parameters). It helps to always test the objective procedure at a single candidate point, as sanity-check. Sorry if I misunderstand you here.

3) There's a battle between the numerical optimizer's wish to establish an optimum that satisfies its optimalitytolerance control option's value versus the numeric integrator's specification to only provide a result that needs to be convergent up to its own accuracy tolerance option epsilon. I know that sounds like a mouthful, but analogous struggles can happen between evalf/Int, fsolve, and Minimize. The "outer" one needs the "inner" one to be accurate enough and stable enough that it can establish its own criteria for convergence. So, in your later example I have forced values for those options. Sometimes it requires fiddling, also because evalf/Int might also require high enough working precision to attain its target accuracy.

Unfortunately, I lost patience fiddling and threw in a scaling factor, to get the integrand to a nice range. It works, but muddies the question of what is really needed.

I have used unapply within the call to evalf(Int(...)) since I've found that all too often it is a useful way to prevent the numeric integrator's zealous approach to checking for discontinuities when presented an integrand in expression form. I didn't check whether it was completly necessary here.

4) I didn't accomodate the case that the numeric integration ever fails to produce a numeric value. The proc3 could (and probably should) test that the numeric integration has returned a numeric value. If it hasn't then it your luck depends on which minimizer you're using, and how it handles ever getting a nonnumeric value for the objective. Some won't handle the discontinuity that a default "large" fallback value can cause. Some won't handle a nonnumeric Float(undefined). DirectSearch handles occasional nonnumeric values from the objective better. 

5) Keep in mind that for the multivariate case Optimization is a local rather than a global optimization package. (The add-on package DirectSearch provides global optimizers.)

The attached worksheet runs pretty quickly. But I had to fiddle to get that.

Download MaplePrimes_Minimization_ac.mw

restart;

expr := y[1](t)*y[2](t)*(y[1](t)+y[2](t))^3:

frontend(diff, [expr, y[1](t)]);

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

Physics:-diff(expr, y[1](t));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

subs(y1t=y[1](t),diff(subs(y[1](t)=y1t,expr),y1t));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

thaw(diff(subs(y[1](t)=freeze(y[1](t)),expr),freeze(y[1](t))));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))