There are quite a few "lesser" problems with the code (loop summation instead of using the add command, inverting the same Matrix many times, repeated computation inside loops of sub-terms which don't depend on the loop-index, etc).

But as far as can see one of the primary efficiency problems is the cost of computing MatrixPower(V, 1/2) for V a float Matrix. Examination indicated to me that was the primary reason why Digits had to be set so high, and also that it was incurring a large port of the time cost.

It turns out that (at present) MatrixPower uses an inaccurate way to compute the square-root of a float Matrix, since (at present) it is just getting redirected to LinearAlgebra:-MatrixFunction to use a general purpose algorithm.

And the only thing that is being done with Vx the computed Matrix-square-root of V is obtaining and utilizing its eigenvalues and eigenvectors. Now, I don't know what your algorithm is supposed to accomplish overall. But for this particular subtask it seems to me that you could do away entirely with that MatrixPower(V, 1/2) call, and compute the eigenvectors and eigenvalues of the Matrix-square-root of V much more directly.

Here's my modification that makes this change. See the comments where the code changed.

FP3_modfidea.mw

This speeds things up enough that the loop-computations have changed from being relatively incidental to now contributing some of the greatest efficiency costs. So I suspect that the overall computation can be made faster still, but  I don't have time to work on it right now.

If there are still places where Digits needs to be raised, for higher N, then certainly that could be done. But such an increase to the working precision should probably only be done for the the duration of the subcomputations where it is necessary. Set Digits back down again, where it doesn't matter.

It's even possible that the eigen-solving might require a higher value for Digits as N gets large. But even then it'll be faster to utilize the simpler Eigenvectors(V) results than to call MatrixPower.

I'm going to log a bug report, that MatrixPower is not doing the best thing itself in this float Matrix example. Note that even if MatrixPower did the right thing here (accurately) it would still be the wrong thing to use: it'd make Maple compute eigenvectors twice, when once is adequate.

acer

I tried to speed it up, but a smoothish curve for u=0.1 .. 10 still took about an hour on my i7.

The plots look a lot smoother when viewed in the GUI. (Why is the plot rendering so poor for uploaded workeets on this site?)

I changed the first (of the innermost pair) integral to go from 1 to infinity (instead of 100), as that allowed exact integration after a changed of variables and assumptions that matched the outer integration variables' ranges. Let me know it this is not ok. (The integrand is decaying...)

It was a bit awkward to make some of the integrands be procedures, but I found that I needed to avoid the costly expansion of integrands that evalf/Int attempts. So instead I used the expandoff command to prevent such unwanted expansions. (Without disabling the relevant expansions the computing of the outer integral caused memory exhaustion on my computer.)

I found the final computed plot to be strange over the 0.0 .. 0.1 domain. The are also some wiggles over the 0.1 .. 2.0 domain. The latter I tried to resolve by raising Digits (15, and 30, as attempts), and refining the various epsilon tolerances, and replacement of the NAG methods by _Dexp. But that just made things slower, without the plot being smoother over than subdomain.

Note that repeating the exact integration result, after a full run through the worksheet, needs either the right call to expandon or a restart.

Download tripplint_expandoffD.mw

acer

Previous answers are fine, though this particular plotting example can also be handled by rescaling the expression. And below is another way to avoid "hardware floats".

Download sfloats.mw

acer

This is what I have so far...

It produces a much simpler form of the result than what you previously obtained from int, and it demonstrates that this is equal to your expected form (involving coth).

It also shows how your short form involving coth can be turned into this result (much simpler than previously obtained) programatically.

But so far I haven't found a way to turn this new (simpler) form of the result into your equivalent form involving coth.

Download trighint1.mw

acer

restart;

X:=Vector([seq(-evalhf(2*Pi)..evalhf(2*Pi),.1)],datatype=float[8]):

Y:=map(sin,X):

F:=proc(x, {degree::posint:=3})
    if not type(x,numeric) then return 'procname'(args); end if;
    CurveFitting:-ArrayInterpolation(X, Y, x,
                     ':-method'=':-spline', ':-degree'=degree):
end proc:

evalf(Int(F, 0..Pi));
                   1.9999997215809067

evalf(Int(x->F(x,degree=3), 0..Pi));
                   1.9999997215809067

evalf(Int(x->F(x,degree=5), 0..Pi));
                   1.999999999903325

acer

With a slightly stronger hint for the restriction on parameter s2 the Optimization:-LSSolve command can also be used here.

Download data_fit_B.mw

acer

You didn't upload your data file (you can just zip it, and fake the extension, btw...) so I have to guess as to how you made and saved it.

But here is the result of a little experiment, using 64bit Maple 17 for Windows 7.

restart;

# Import it, with all entries imported as strings.
#
M := ImportMatrix(cat(kernelopts(homedir),"/My Documents/mapleprimes/md1.csv"),
                  delimiter=",",datatype=string);

                      ["Symbol"  "2015-08-15"  "2016-08-15"]
                      [                                    ]
                 M := [ "TFSC"      "9.76"        "9.62"   ]
                      [                                    ]
                      [ "PIH"       "7.63"        "7.57"   ]

# Strip off the restrictive string datatype.
#
M := Matrix(M, datatype=anything):

# Parse everything but the first column and row, to get numbers.
#
M[2..,2..]:=map(parse,M[2..,2..]):

M;

                   ["Symbol"  "2015-08-15"  "2016-08-15"]
                   [                                    ]
                   [ "TFSC"       9.76          9.62    ]
                   [                                    ]
                   [ "PIH"        7.63          7.57    ]

The file md1.zip contains the .csv file I used.

md1.zip

md1.mw

[edit] Now that the OP has uploaded a zipped .csv file, I checked and the above works the same on that data1.csv file.

acer

@jaypantone Do you see the same buffered/delayed effect for userinfo statements? (Those have the nice feature that you can turn them all on or off with a single command, and also control different levels of detail. They ought to display in a visually similar way in CLI and Std GUI)

Eg,

p := proc()
  userinfo( 1, p, `This is a userinfo statement`);
  userinfo( 3, p, `This is a level 3 a userinfo statement`);
  userinfo( 1, p, nprintf("This is a level %a userinfo statement,\nconstructed with nprintf", 1));
end proc:

p();

infolevel[p]:=3:
p();

infolevel[p]:=1:
p();

It's interesting that you don't see the delay problem with the `print` command. Here's a thought: does using `print` or `lprint` to print a single character string such as " ", interleaved with your delay-buffered `printf` statements, help at all?

Is this the kind of thing you want?

restart:

x:=[1992,1994,1996,1998,2000,2002,2004]:
y:=[325,310,266,217,178,159,161]:

Statistics[ColumnGraph](y,captions=[seq(sprintf("%a\n\n",_n),_n=x)]);

kernelopts(version);

          Maple 12.02, X86 64 WINDOWS, Dec 10 2008 Build ID 377066

acer

f := 5*A^4*B^3*C^7   +   3*A^2*B^1*C^7  +  31*A^3*B^6*C^11;

                          3  6  11      4  3  7      2    7
                 f := 31 A  B  C   + 5 A  B  C  + 3 A  B C 

Exponents := proc(p,x)
  local F;
  if p::`+` then
    [seq(degree(F,x), F in p)];
  else
    [degree(p,x)];
  end if;
end proc:

Exponents(f, A);
                                  [3, 4, 2]

Exponents(f, B);
                                  [6, 3, 1]

Exponents(f, C);
                                 [11, 7, 7]

Exponents(31*A^3*B^6*C^11, C);
                                    [11]

acer

If you want to call the LinearAlgebra package's commands with their short form names then at the start of your procedures include the statement,

uses LinearAlgebra;

For example,

f := proc(A)
  local res;
  uses LinearAlgebra;
  res := Basis(A);
end proc:

f( [<0,1,1>,<1,2,1>,<1,1,0>] );

                                 [[0]  [1]]
                                 [[ ]  [ ]]
                                 [[1], [2]]
                                 [[ ]  [ ]]
                                 [[1]  [1]]

If you are writing a module based package then you can do that within each of its procedures, be they exports or local (to the module) procedures.

Basically the `uses` statement sets up the relevant name bindings within your procedure, somewhat analogously to how `with` works at the top-level. (nb. `with` is not the thing to do, inside a procedure.)

The `uses` facility is documented in the Description section of the help-page for proc (procedure).

acer

restart:

limit((b^(d+1)-1)/((b-1)*b^d), d = infinity) assuming b>1;

                                      b  
                                    -----
                                    b - 1

limit((b^(d+1)-1)/((b-1)*b^d), d = infinity) assuming b<-1;  # hmm

                                  undefined

Q := simplify( (b^(d+1)-1)/((b-1)*b^d) ) assuming b<>1;

                                         (-d)
                                    b - b    
                               Q := ---------
                                      b - 1  

R := limit(Q, d = infinity) assuming b<-1 or b>1;

                                        b  
                                 R := -----
                                      b - 1

limit(Q, d = infinity) assuming b>=0 and b<1;

                                  infinity

limit(Q, d = infinity) assuming b>=-1 and b<0;

                                  undefined

What is going on as b->1 and d->infinity?

acer

Do you mean something like this? I am guessing that you want subscripts rather than underscores, which in your Maple 2015 can be done using a literal double-underscore.

plot( [[0,0]], labels=[Omega/omega__n,a__0(m)] );

If you want that m to be upright and not italic, you could do it as say,

labels=[typeset('Omega/omega__n'), typeset('a__0(`#mn("m")`)')]

acer

Is this the kind of thing you mean?

restart:

N := 50:

e2 := [seq(k*Pi/N, k=-N..N)]:
e3 := evalf(map(sin, e2)):

e4 := CurveFitting:-Spline(e2, e3, x):

PD := proc(p, v)
  local L;
  L:=convert(p, list);
  piecewise(seq([L[2*i-1],diff(L[2*i],v)][], i=1..floor(nops(L)/2)),
                diff(L[-1],v));
end proc:

de4 := PD(e4, x):

plot([e4, de4], x=-Pi..Pi);

acer

There is a call in the OpenMaple API for querying the interrupt signal. A lot of dsolve/numeric is in external compiled code, which could make use of that.

acer