acer

20490 Reputation

29 Badges

15 years, 30 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Carl Love 

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.

@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.

[note. I have removed the link that I added. The OP can add a correct link and attribution if he wishes.]

 

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

@Josolumoh 

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.

Where is the code itself?

@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 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...).

restart;

kernelopts(version);

`Maple 2018.1, X86 64 LINUX, Jun 8 2018, Build ID 1321769`

with(LinearAlgebra):

M:=parse(FileTools[Text][ReadFile]("1.txt")):

Dimensions( M );

70, 70

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]);;

Vector(5, {(1) = -73328119.29, (2) = -19406541.76, (3) = -216657.530800000008, (4) = -12263.6964700000008, (5) = -975.701514900000006}), Vector(5, {(1) = -73328119.29, (2) = -19406541.76, (3) = -216657.5308, (4) = -12263.69647, (5) = -975.7015149})

infolevel[LinearAlgebra]:=0:

Norm( Epd - sort(Re(E)) );

0.121761997198733919337099580190866e-107

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

48, 0.30e-5

49, 0.18e-6

50, 0.13e-7

51, 0.14e-8

52, 0.16e-9

53, 0.18e-10

54, 0.14e-11

55, 0.16e-12

56, 0.25e-13

57, 0.17e-14

58, 0.11e-16

 

Download evals_gen_sym_posdef.mw

@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)

Put your followup queries on this topic as comments here in this same thread, please.

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.)

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

@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. 

@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));

[1., 6.573224507, 4.850242691, 3.616575701, 2.840981033, 2.493303585, 2.424667056, 2.422212016, 2.422208959, 2.422208959, 2.422208959]

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

[-8., -8.487501254, -8.568487419, -8.570502786, -8.570504010, -8.570504009, -8.570504009, -8.570504009, -8.570504009, -8.570504009, -8.570504009]

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

S0*k^2*A(y)^(10/3)/(Q*abs(Q)*P(y)^(4/3))-1

dee := diff(ee,y);

-(4/3)*S0*k^2*A(y)^(10/3)*(diff(P(y), y))/(Q*abs(Q)*P(y)^(7/3))+(10/3)*S0*k^2*A(y)^(7/3)*(diff(A(y), y))/(Q*abs(Q)*P(y)^(4/3))

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

(2/3)*S0*k^2*A(y)^(10/3)*(5*(diff(A(y), y))/A(y)-2*(diff(P(y), y))/P(y))/(Q*abs(Q)*P(y)^(4/3))

simplify( hdee -dee );

0

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]);

(9/25000)*((y+12)*y)^(10/3)/(12+2*y*2^(1/2))^(4/3)-1

Digits := 30;

30

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

2.42220895913004814986049595028

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

-0.4597438638e-29

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

-8.57050400897802308341902137162

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

-0.3528983078e-29-0.9378067234e-100*I

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

-(9/100000)*y^3*(y+12)^3*(-(4*y+48)*y*(6+y*2^(1/2))^2)^(1/3)/(6+y*2^(1/2))^2-1

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

-8.57050400897802308341902137162

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

S0*k^2*((m*y+b)*y)^(10/3)/(Q*abs(Q)*(b+2*y*(m^2+1)^(1/2))^(4/3))-1

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

S0*k^2*(-m*y^2-b*y)^(10/3)/(Q*abs(Q)*(-b-2*y*(m^2+1)^(1/2))^(4/3))-1

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

-8.57050400897802308341902137162

 

Download newtthing.mw

@MAJP Each iteration through the loop , do you want to be defining y[i] in term of the previous value y[i-1] ?

Trying to assign y[i] in terms of a formula with y[i] is not right.

You'll  also need to consider whether the other indexed names in the formulas in the loop should use the current iteration's values indexed by i, or the previous iteration's values indexed by i-1.

I feel sure that, by thinking about this aspect, you'll understand what the iteration is supposed to accomplish.

 

@sand15 

Whenever I install Maple the next thing I do is edit the properties of the launcher, to add (along with a blank space as separator) the option
   -standalone
to the end of the invocation of the GUI executable.

This has the effect that I can launch separate worksheets under separate Java jvm processes. And if one of them stalls or freezes or crashes (in the GUI sense) then the other one(s) stay ok.

It'll even let me open the same worksheet twice.

Memory is cheap. Back around Maple 10 the few hundred MB footprint of the GUI+jvm seemed huge. Now it seems modest.
 

First 72 73 74 75 76 77 78 Last Page 74 of 439