## 20821 Reputation

15 years, 57 days

## trig...

@Rim So, clearly you need at least some amount of trig simplification, because your examples might include the need to make, say, sin(t)^2+cos(t)^2 -> 1.

If you're building up bigger and bigger expressions within a loop, it's might be best to keep the intermediate expressions "simplified" at each iteration. But you may wish to avoid everything except basic factoring/grouping and that trig simplification. That is, you may be better off with some degree of control about expression swell and trig simplification at each iteration, but without `simplify` trying its hardest to get the very "simplest" representation at each iteration.

I'm mostly just guessing, because I don't have the full package or the complete examples.

So I'll make a few suggestions for replacement to the costly (according to your analysis) call to just `simplify`. You could try them out separately and re-profile. Perhaps one of them might be faster for your full examples.

The idea is to replace the call  simplfy(ee)  which you've identified as a bottleneck, where ee is just the single expression argument.

simplify(ee, trig)
simplify(normal(ee, trig))
`simplify/nosize`(ee, trig)
`simplify/nosize`(normal(ee),trig)

There are other actions that might be thrown into this mix (combine, factor, frontend, etc).

What might be even more fruitful is if you could produce some enormous, "wrong" final result that has not been simplified as you need. If you can produce such, you could upload a worksheet here as an attachment (big green up-arrow).

What version are you using, BTW?

some_trig_simp.mw

## suggestion...

@Rim Suppose that the code just calls simplify at some internal part of the code.

One possibility is that the code does this inefficiently, eg. calling it unnecessarily and repetitively inside a loop. Presumably you are already checking out for that.

Another possibility is that simplify is too costly a tool for the simplification task at hand. We might be able to help there, but it would likely require several large examples of both input expression and desired output expression. There may be a combination of lighter weight commands that produce the same simplifications.

If you repost the procedure's source code here the please provide either a URL where you got it or the name of the author(s).

## published where?...

Free-dealing/Free-use exceptions to copyright for some purposes can include the need for proper attribution.

But also (and if applicable here then supercedent) copyright free-use pertains to published work. If the code itself has not been published by its rightful owner then the question of copying/reproducing it is moot. Publishing someone else unpublished work is not copying and is not ok on the fair-dealing/fair-use grounds governing copyrights.

For example, this thesis says, in its Appendix A, that some of the "files can be received from the Chair for Automation and Control of the Department for Measurement and Automation of the University of the German Armed Forces in Munich". That need not mean that the code available only in said files is published.

## someone else's work...

@Carl Love He previously posted the full code of the procedure skewSimplifier. I deleted it as an unattributed procedure from someone's else's package and thesis, posted without the author's explicit permission.

If he can show the permission of the author of that full procedure to re-post it here, or at the very least add proper attribution, then I'll let it stand.

Even fair use (or fair-dealing) exceptions to copyright generally requires proper attribution of authorship.

## why Maplets?...

I'm interested to know why you want to program this using Maplets rather than Embedded Components. (Each has its merits, and drawbacks.)

## Matrix...

These both work for me as I might expect,

```ExportMatrix("EM.txt", Matrix(Vector[column](5,i->i)), target = csv);

ExportMatrix("EM.txt", Matrix(Vector[row](5,i->i)), target = csv);
```

So perhaps you could try wrapping your result (from Statistics:-Sample) in a call to Matrix. (A column Vector or a row Vector produces the nx1 or 1xn Matrix, accordingly.) This seems to work reasonably with ExportMatrix.

## the code...

Where is the code itself?

## clues...

@tomleslie The preamble in the OP's code, looked to me as if it were being used in the construction of the eqs[i], but unfortunately not in a fully programmatic way (or, at least, not a way revealed to us). The preamble code made me suspicious that the equations might be underdetermined and possibly admit infinitely many solutions. Impossible to tell for sure without getting rid of the floats and discovering just how the eqs[i] are built. But I'd already noticed several familiar floats.

It's a real shame that so many posters here only reveal a portion of their goal. It is a very common phenomenon that a Questioner here will only ask about the problematic methodology that they've chosen for just one small step in a (possibly dubious) overall approach to a wider problem.

## why floats?...

Why do you use floating-point numbers in the first place?!

I mean, are you trying to use these values, or something similar?

```phi[1][0], phi[2][0], phi[1][1], phi[2][1], phi[1][2], phi[2][2]
:= sqrt(2), sqrt(2), (4*t-1)*sqrt(6), (4*t-3)*sqrt(6),
(1/2)*sqrt(10)*(3*(4*t-1)^2-1), (1/2)*sqrt(10)*(3*(4*t-3)^3-1)```

If so then why not use exact rationals for alpha, lambda[3], etc. And why not subsequently construct the eqs[i] from the X.Y,H[1],etc and their derivatives?

Maybe you'd be able to get an exact root, which you could approximate quickly and with less accumulated roundoff error.

Maybe you'd even be able to get an exact symbolic formula for a generic form of the root for general k and M. (I can't tell, because you've already done the dubious step of cut and pasting just to get the formulas used used to construct eq[i]. And the formulation has been obscured...).

## A symmetric, B symmetric positive defini...

 > restart;
 > kernelopts(version);

 > with(LinearAlgebra):
 > Dimensions( M );

 > interface(rtablesize=70):
 > Digits:=150:
 > infolevel[LinearAlgebra]:=1:
 > A:=Matrix(eval(M,P=0),shape=symmetric):
 > B:=Matrix(-diff(M,P), shape=symmetric):
 > E:=CodeTools:-Usage( Eigenvalues(A,B) ):

Eigenvalues: calling external function

Eigenvalues: CLAPACK sw_dggevx_

memory used=4.10GiB, alloc change=38.00MiB, cpu time=15.01s, real time=12.95s, gc time=4.82s

 > Bpd:=Matrix(B,shape=symmetric,attributes=[positive_definite]):
 > Epd := CodeTools:-Usage( Eigenvalues(A,Bpd) ):

Eigenvalues: calling external function

Eigenvalues: CLAPACK sw_dsygvd_

memory used=498.11MiB, alloc change=-4.00MiB, cpu time=1.65s, real time=1.39s, gc time=567.68ms

 > evalf[10](Epd[1..5]), evalf[10](sort(Re(E))[1..5]);;

 > infolevel[LinearAlgebra]:=0:
 > Norm( Epd - sort(Re(E)) );

 > for d from 48 to 58 do   Digits:=d:   Epd2 := Eigenvalues(A,Bpd):   nrm := Norm(Epd - Epd2):   print(d, evalf[2](nrm)); end do:

 >

## mechanism is weak...

@rquirt Yes, thanks for the upload. (You don't need to insert all the content of the upload, every time. Inserting the link suffices.) But thanks.

The mechanism the fsolve uses to determine what are the units (of both the input and the output) is not strong. It seems to require the literal presence of the output units in (some addend of) the expression to be solved. But that is quite fragile and breaks for some choices of fsolve calling sequence. It is inadequate for both operator form as well as the unevaluated function call. I will submit some bug reports.

Here are two workaround that succeed for your procedure test.

The first adds an insignificantly small amount of Unit(W), which is adequate to let fsolve figure out the output unit and not balk later. I note that in your examples with just the single function call you already had a constant term in Unit(W). It was the absence of such which "broke" the test example.

The second, which was also in my answer, essentially strips out input and output units and proceeds with raw numbers.

```fsolve(test(x)+0.1e-99*Unit(W), x = 20*Unit(degC) .. 30*Unit(degC));

23.94824265 Unit(°C)

fsolve(x -> test(x*Unit(degC))/Unit(W), 20 .. 30)*Unit(degC);

23.94824265 Unit(°C)
```

## followups...

Please upload .mw worksheets (green up arrow in the Mapleprimes editor) since you're having copy and paste issues. (And never paste in 2D input code as plaintext here.)

## Actually...

@Christopher2222 Actually you should describe your actual background problem up front.

## several changes...

@rquirt Units:-Simple allows you to add blah*Unit(W) and the unevaluated function call to your proc (which itself returns a quantity with units).

I removed the call to `with(Units:-Standard)` at the top level. The procedure revision i gave works without it. I utilized `uses` in the proc instead. (It is not great programming to write a proc that only works if `with`is called outside it.)

I also made your proc return unevaluated when the first argument is not numeric*unit or numeric.

There was also a left- right- quote match typo in your original.

Better to upload actual working code or worksheet.

## typo...

@MAJP Firstly, I commend you for putting in real effort at solving this. It's refreshing to see.

I suspect that you have made two or three typo's in writing down the formula for Fdiff[i].

The first is that you are missing a multiplication sign between two pairs of round brackets. So that ends up being parsed as function application rather than multiplication. (Maple allows you to apply a floating-point number to something else, and it returns that same float.)

The second seems to be that you have (F[i]-1) instead of (F[i]+1) in the formula for Fdiff[i].

Also, as Carl mentioned, you have an y[1] which ought to be y[i] in the formula for dAdy[i].

If I make those adjustments then the derivative of F appears to match Fdiff. And then the starting point y0=1 allows convergence to the root near y=2.4222.

In order to get convergence to the apparent root near y=-8.57 I wrapped the formula for y[i+1] with two commands whose purpose is to remove any small imaginary component (where small means close to the working precision.) You can get rid of that if you find it's out of scope for your task.

I suppose you realize that you don't really need to keep track of all the values of A,P,Fdiff,etc for all i values. You could just as well use non-indexed names for those. The indexing of Y is key to your approach, of course.

 > restart;
 > Newton := proc (n, y0, Q, b, m, k, S0)   local A, P, Sf, y, F, dAdy, dPdy, Fdiff, i, ylist;   y[0] := evalf(y0);   for i from 0 to n do     A[i] := (b+m*y[i])*y[i];     P[i] := evalf(b+2*y[i]*sqrt(1+m^2));     Sf[i] := evalf(Q*abs(Q)*P[i]^(4/3)/(k^2*A[i]^(10/3)));     F[i] := S0/Sf[i]-1;     dAdy[i] := b+2*m*y[i];     dPdy[i] := evalf(2*sqrt(1+m^2));     #Fdiff[i] := -(4/3)*S0*k^2*A[i]^(10/3)*(dPdy[i])/(Q*abs(Q)*P[i]^(7/3))     #            +(10/3)*S0*k^2*A[i]^(7/3)*(dAdy[i])/(Q*abs(Q)*P[i]^(4/3));     Fdiff[i] := (2/3)*(F[i]+1)*(5*dAdy[i]/A[i]-2*dPdy[i]/P[i]);     y[i+1] := simplify( fnormal( y[i]-F[i]/Fdiff[i] ), ':-zero' );   end do;   ylist := [seq(y[i], i = 0 .. n)];   return ylist end proc:
 >
 > Newton(10, 1, 5, 12, 1, 30, 10^(-5));

 > Newton(10, -8, 5, 12, 1, 30, 10^(-5));

 > ee := S0/(Q*abs(Q)*P(y)^(4/3)/(k^2*A(y)^(10/3)))-1;

 > dee := diff(ee,y);

 > hdee := (2/3)*(ee+1)*(5*diff(A(y),y)/A(y)-2*diff(P(y),y)/P(y));

 > simplify( hdee -dee );

 > this := eval(eval(ee, [A(y)=(b+m*y)*y,                   P(y)=b+2*y*sqrt(1+m^2)]),              [Q=5, k=30, S0=10^(-5), m=1, b=12]);

 > Digits := 30;

 > sol1 := fsolve(this, y=1..10);

 > evalf[100](eval(this, y=sol1)): evalf[10](%);

 > sol2 := fsolve(abs(this), y=-9..-8);

 > evalf[100](eval(this, y=sol2)): evalf[10](%);

 > foo := simplify(combine(expand(this)),size) assuming y>-9, y<-8;

 > fsolve(foo, y=-9..-8);

 > bar := eval(ee, [A(y)=(b+m*y)*y,                   P(y)=b+2*y*sqrt(1+m^2)]);

 > bbar := simplify(evalc(bar),size) assuming real, m=10, b<13, y>-9, y<-7;

 > fsolve(eval(bbar, [Q=5, k=30, S0=10^(-5), m=1, b=12]), y=-9..-8);

 >