acer

33400 Reputation

29 Badges

20 years, 331 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

@Carl Love Yes, thanks. There is also a pair DDTTRFB/DDTTRSB in MKL for factoring and backsolving for diagonally dominant tridiagonal matrix. (That factorization returns no pivot list, only using dl,d,du. I do not know what it uses, eg. Thomas' algorithm.)

Like you I consider it worthwhile to be able to separate the factorization from the solving step (for a given RHS). For the LU with pivoting that might be possible as a special list return from LUDecomposition, from the other tridiagonal solver I described in my Answer. That prefactorization list might then be passed into LinearSolve. That's how it works now for Cholesky, QR, and LU (but not band LU). It's a little awkward, and all those routines have overhead as they figure out how to dispatch, whether it's a float or exact problem, etc.

As an alternative... perhaps it would be worthwhile to create BLAS and LAPACK packages, with all the in-place and other goodness. They could run hardware and software floats, as possible. But a goal would be low overhead. It would also offer some nice efficiencies like being able to specify by semaphore that a Matrix argument was to be taken as "transposed", without copying or computation cost, eg. solve A^%T . x = b, or multiply a*A^%T . B , etc. I wouldn't be surprised if much of the code could be auto-generated. Another benefit is that CodeGeneration[C] might be taught about it. For some time I've thought about this.

@Rouben Rostamian  Thanks.

And yes, indeed, for a different example I might have had to call frontend differently (which is a common phenomenon for many frontend uses...).

So Physics:-diff is a more general solution for this kind of task.

note. It's often a good idea to allow these (default expression types, {`+`,`*`}) through as unfrozen,

  frontend(diff, [x(t)^2 + sin(x(t)+y(t)), x(t)], [{`+`,`*`},{sin}]);
                            2 x(t) + cos(x(t) + y(t))

  frontend(diff, [x(t)^2 + sin(x(t)+y(t)), x(t)], [{`+`,`*`,specfunc(sin)},{}]);
                            2 x(t) + cos(x(t) + y(t))
As the examples gets more involved and confusing even targeted freeze/thaw starts to look somewhat attractive. Having to double-check every answer robs an approach of merit. I guess I ought to have mentioned Physics:-diff first.

@Carl Love This member usually posts worksheets using Maple 18, and that might still be his working version. Your cited post's code would need revising, to run in that version, I believe.

Fyi, when I ramped up the size and number of RHSs to both be very large I sometimes obtained a result that contained Float(undefined).  However I was sometimes able to obtain a solution by other means in that case, ie. I'm not sure that I was encountering inconsistency in the random example. It was rare. I didn't poke around.

Timing performance of your prior code seemed quite similar to that of the LAPACK tridiagonal solver (as I accessed it), although the latter usually produced a slightly smaller residual. The timing differences were affected by the order in which they were run, and aspects that hinted the loading of shared external objects was a factor.

@Carl Love Yes, sorry.

@Carl Love Right you are, of course (I knew that but was not thinking, it might be a WFH thing).

I think that it's interesting how much garbage collection is part of some of these timings.

Something that I was considering mentioning was that these kinds of schemes can vary in performance according to which aspect of the problems grows in size. So for this problem it might be whether the number of items or the number of keys (or both) grows large.

When I was writing the simple loop code I showed in my Answer I was primarily thinking about dealing with a small number of unsorted keys relative to the number of items. I would write it differently for a large number of keys and, as I did mention in my Answer, I'd consider using ListTools:-Classify. (I lnow I used an example with 10^3 keys in comparison with Kitonum's code -- what I do and what I think don't always match.)

For a large enough number of items but a small number of keys (here, distinct eigenvalues lists) the scheme I coded performs modestly better than Carl's scheme. 

But Carl's scheme scales up much better as both the number of keys as well as the number of items grows large.

Another consideration which sometimes comes into play is how many such computations need to be done, and whether any part is repeated.

And that is an important takeaway, I think -- that the scheme should be taken which best suits the character of the actual data.

Unfortunately, in this question we have not been provided with the actual data, so there is no way to know what would be best except for guessing.

We can also observe that, for some sizes of problem, the repetition and ordering of the compared computations matters and has an effect. Therefore a better comparison would involve separate sessions. (The generated examples below are random and hence session-dependent, but should be large enough that the timings are representative and still meaningfully comparable across sessions.)

Here are some more data points, for fun:

restart;

#Acer's code (verbatim):
#
GL := proc(LL)
  local U,C,T2,t,u,i,L;
  T2 := {map[2](op,2,LL)[]};
  for i from 1 to nops(T2) do
    C[T2[i]] := 0:
  end do:
  for L in LL do
    for t in T2 do
      if L[2]=t then
        C[t] := C[t] + 1;
        U[t][C[t]] := L[1];
        break;
      end if;
    end do;
  end do;
  [seq([lhs(u),convert(rhs(u),list)],u=op([1,..],U))];
end proc:
#

N1,N2:= 10^1,10^6:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GL(MyNestedList)):

memory used=459.37MiB, alloc change=106.02MiB, cpu time=11.55s, real time=6.49s, gc time=7.11s

N1,N2:= 5,10^6:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GL(MyNestedList)):

memory used=398.44MiB, alloc change=22.91MiB, cpu time=7.52s, real time=5.12s, gc time=3.37s

N1,N2:= 10^2,10^5:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GL(MyNestedList)):

memory used=151.74MiB, alloc change=-45.81MiB, cpu time=2.60s, real time=1.83s, gc time=1.13s

N1,N2:= 10^4,10^4:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GL(MyNestedList)):

memory used=0.73GiB, alloc change=-7.64MiB, cpu time=15.83s, real time=13.14s, gc time=4.03s

 

 

restart;

#Carl's alternative:
#
GraphEigenvalueClassify:= (L::list([{Graph, name}, list]))->
    [lhs,rhs]~(
        (op@map)(
            g-> [map(`?[]`, g, [1])[]],
            ListTools:-Classify(g-> g[2], L)
        )
    )
:

N1,N2:= 10^1,10^6:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GraphEigenvalueClassify(MyNestedList)):

memory used=0.94GiB, alloc change=179.12MiB, cpu time=28.27s, real time=16.38s, gc time=18.10s

N1,N2:= 5,10^6:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GraphEigenvalueClassify(MyNestedList)):

memory used=0.52GiB, alloc change=47.28MiB, cpu time=9.87s, real time=7.67s, gc time=3.40s

N1,N2:= 10^2,10^5:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GraphEigenvalueClassify(MyNestedList)):

memory used=55.10MiB, alloc change=0 bytes, cpu time=576.00ms, real time=576.00ms, gc time=0ns

N1,N2:= 10^4,10^4:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

CodeTools:-Usage(GraphEigenvalueClassify(MyNestedList)):

memory used=21.54MiB, alloc change=0 bytes, cpu time=95.00ms, real time=95.00ms, gc time=0ns

 

Download grouplistcomp.mw

ps. I know that the code snippet to generate random examples could be improved. I didn't consider it important enough to optimize.

@hakan So, do you mean like this?

I'll also put in Kitonum's procedure, in case you wish to compare. Of course the order of the sublists may differ in the result, but it's not difficult to use a list rather than a set of the keys -- arranged as you prefer.

I repeat each twice for the larger example, interleaved, to show that the timing's are not just a consequence of which runs first.

restart;

MyNestedList := [ [g1, [1,2,3]], [g2, [0,1,3]], [g3, [1,2,3]], [g4, [0,1,3]],
                  [g5, [2,4,7]], [g6, [0,1,3]], [g7, [2,4,7]], [g8, [0,-1,9]] ]:

 

GL := proc(LL)
  local U,C,T2,t,u,i,L;
  T2 := {map[2](op,2,LL)[]};
  for i from 1 to nops(T2) do
    C[T2[i]] := 0:
  end do:
  for L in LL do
    for t in T2 do
      if L[2]=t then
        C[t] := C[t] + 1;
        U[t][C[t]] := L[1];
        break;
      end if;
    end do;
  end do;
  [seq([lhs(u),convert(rhs(u),list)],u=op([1,..],U))];
end proc:

 

Ans1 := CodeTools:-Usage( GL(MyNestedList) );

memory used=33.94KiB, alloc change=0 bytes, cpu time=1000.00us, real time=0ns, gc time=0ns

[[[2, 4, 7], [g5, g7]], [[0, -1, 9], [g8]], [[1, 2, 3], [g1, g3]], [[0, 1, 3], [g2, g4, g6]]]

Grouping:=proc(L::listlist)
local L1;
uses ListTools;
L1:=[Categorize((x,y)->x[2]=y[2], L)];
map(p->[p[1,2],map(l->l[1],p)], L1);
#map(p->[p[1,2],map[2](op,1,p)], L1);
end proc:

 

Ans2 := CodeTools:-Usage( Grouping(MyNestedList) );

memory used=0.75MiB, alloc change=0 bytes, cpu time=11.00ms, real time=11.00ms, gc time=0ns

[[[1, 2, 3], [g1, g3]], [[0, 1, 3], [g2, g4, g6]], [[2, 4, 7], [g5, g7]], [[0, -1, 9], [g8]]]

 

nops(Ans1), nops(Ans2), evalb( {Ans1[]} = {Ans2[]} );

4, 4, true

 

N1,N2 := 10^3,10^4:
M := [seq(convert(LinearAlgebra:-RandomVector(3),list),i=1..N1)]:
f:=rand(1..N1):
MyNestedList := [seq([g||i,M[f()]],i=1..N2)]:

 

CodeTools:-Usage( GL(MyNestedList) ):

memory used=135.80MiB, alloc change=48.00MiB, cpu time=2.08s, real time=1.92s, gc time=395.52ms

CodeTools:-Usage( Grouping(MyNestedList) ):

memory used=0.58GiB, alloc change=-4.00MiB, cpu time=8.75s, real time=8.02s, gc time=1.53s

CodeTools:-Usage( GL(MyNestedList) ):

memory used=119.80MiB, alloc change=0 bytes, cpu time=1.72s, real time=1.59s, gc time=275.15ms

CodeTools:-Usage( Grouping(MyNestedList) ):

memory used=0.58GiB, alloc change=0 bytes, cpu time=8.41s, real time=7.66s, gc time=1.57s

 

Download grouplist.mw

@fmoraga So did you not read my earlier response about syntax for multiplication? Making the simple edit of adding explicition multiplication symbols in two places (Eq2, Eq3) seems to clear up the problem.

It is unclear why you prefer one representation of W_dot over others. Perhaps you have reasons for that.

If you want also to eliminate {Q_dot_C,Q_dot_H,T_c,T_h} specifically alongside {W_dot} then you could call eliminate against the union of those together -- and obtain no restrictions on the remaining variables. Or you could simply call solve directly against that union -- since the remaining variables are all free.

And you could then use eval as the second step to extract only the W_dot formula, which is simpler than having another solve call as the second step.

I suppose that you realize that you could also eliminate {Q_dot_C, Q_dot_H, eta} specifically in addition to W_dot. Or, you could simply take Eq1 as the answer, or use only Eq4.

Here are your original elimination, and a couple of others.

restart

Eq1 := W_dot = Q_dot_H-Q_dot_C;

W_dot = Q_dot_H-Q_dot_C

Eq2 := Q_dot_H = UA_H*(T_H-T_h);

Q_dot_H = UA_H*(T_H-T_h)

Eq3 := Q_dot_C = UA_C*(T_c-T_C);

Q_dot_C = UA_C*(T_c-T_C)

Eq4 := eta = W_dot/Q_dot_H;

eta = W_dot/Q_dot_H

Eq5 := eta = 1-T_c/T_h;

eta = 1-T_c/T_h

List := eliminate({Eq1, Eq2, Eq3, Eq4, Eq5}, {Q_dot_C, Q_dot_H, T_c, T_h}):

solve(List[2], W_dot);

{W_dot = UA_C*UA_H*eta*(T_H*eta+T_C-T_H)/(UA_C*eta+UA_H*eta-UA_C-UA_H)}

another := eliminate({Eq1, Eq2, Eq3, Eq4, Eq5}, {Q_dot_C, Q_dot_H, T_c, T_h, W_dot}):

{}

W_dot = eval(W_dot, another[1]);

W_dot = UA_C*UA_H*eta*(T_H*eta+T_C-T_H)/(UA_C*eta+UA_H*eta-UA_C-UA_H)

anotherS := solve({Eq1, Eq2, Eq3, Eq4, Eq5}, {Q_dot_C, Q_dot_H, T_c, T_h, W_dot}):

W_dot = eval(W_dot, anotherS);

W_dot = UA_C*UA_H*eta*(T_H*eta+T_C-T_H)/(UA_C*eta+UA_H*eta-UA_C-UA_H)

``

alt := eliminate({Eq1, Eq2, Eq3, Eq4, Eq5}, {Q_dot_C, Q_dot_H, T_c, W_dot, eta}):

{}

W_dot = simplify(eval(W_dot, alt[1]), size);

W_dot = (T_H-T_h)*UA_H*(UA_H*(T_H-T_h)+UA_C*(T_C-T_h))/(UA_H*(T_H-T_h)-UA_C*T_h)

new := solve({Eq1, Eq2, Eq3, Eq4, Eq5}, {Q_dot_C, Q_dot_H, T_c, W_dot, eta}):

W_dot = (T_H-T_h)*UA_H*(UA_H*(T_H-T_h)+UA_C*(T_C-T_h))/(UA_H*(T_H-T_h)-UA_C*T_h)

 

Download Tarea1_ac3.mw

@nikoo This answer is nonsense.

Start off by stating clearly what it means to you for a point to be a solution to this system.

Does it mean that the forward error (substituting results into the rhs-lhs of the equations, and taking absolute values, say) must be less than some specific value, and if so what? At what working precision ought that to hold. Be specific.

Are you looking for a value such that -- at high enough working precision, and if the point is refined to enough accurate decimal places -- that forward error could be made arbitrarily small?

Or would you be content with results in which the points only minimize the forward error (or some other suitable objective) and that it can be approached as an optimization problem? You may be forced into this position due to the fact that your equations already have floating-point coefficients(!). Actual "roots" may not exist, because of imprecise coefficients in the equations.

The nature of your final results (approximates actual roots, only minimizes some objective, etc) may well bear on the validity of use of the results in further scientific computation. Few people consider that properly.

@Jameel123 If I go here then I see a link to here, which contains a link to this, which in turn has two links, apparently for version 1 and 2 of a package.

The second link has a download link (a compressed tar file named adyz_v2_0.tar.gz, apparently containing a "DESOLVII" instructional worksheet and a Desolv-IIa.mpl file).

@Lali_miani If you are still using Maple 12 (like many of your previous Questions) then you can do it as follows.

(The elementwise operator `~` did not work in Maple 12. But it is not necessary to use that -- or even to map -- since expressions are taken by solve as equal to zero by default.)

restart;
kernelopts(version);

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

P:=x->x^4+x^3+a*x^2+sqrt(2)*x+b:

evalc([Re,Im](P(1+I)));

                   1/2           1/2
            [-6 + 2    + b, 2 + 2    + 2 a]

solve(%);

                        1/2
                       2              1/2
             {a = -1 - ----, b = 6 - 2   }
                        2

eval(P(x),%);

                /      1/2\
       4    3   |     2   |  2    1/2          1/2
      x  + x  + |-1 - ----| x  + 2    x + 6 - 2
                \      2  /

solve(%);

                                     1/2                        1/2
 1 + I, 1 - I, -3/2 - 1/2 I + 1/2 I 2   , -3/2 + 1/2 I - 1/2 I 2

There is no attachment.

@John Blacksad While these are even more ways (including using commands explictly), you could also get the characteristic matrix or the eigenvalues more directly. (I don't know whether you were doing an assignment question that required you do each step.)

restart

with(LinearAlgebra)

tau := Matrix(3, 3, {(1, 1) = 200, (1, 2) = 0, (1, 3) = 0, (2, 1) = -50, (2, 2) = -800, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 200})"(->)"lambda^3+400*lambda^2-280000*lambda+32000000"(->)"[[lambda = -800], [lambda = 200], [lambda = 200]]

 

tau := Matrix(3, 3, {(1, 1) = 200, (1, 2) = 0, (1, 3) = 0, (2, 1) = -50, (2, 2) = -800, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 200})"(->)"Vector(3, {(1) = -800, (2) = 200, (3) = 200})"(->)"sigma

sigma[1]

-800

sigma[2]

200

sigma[3]

200

``

Download maple_ac2.mw

And here are some ways, using commands explicitly,

restart

with(LinearAlgebra)

``

tau := Matrix(3, 3, {(1, 1) = 200, (1, 2) = 0, (1, 3) = 0, (2, 1) = -50, (2, 2) = -800, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 200})

Matrix(%id = 18446883827713064582)

sigma := Eigenvalues(tau)

Vector[column](%id = 18446883827713055910)

sigma[1]

-800

sigma[2]

200

sigma[3]

200

CharacteristicMatrix(tau, lambda)

Matrix(%id = 18446883827713051934)

Determinant(%)

(lambda-200)^2*(lambda+800)

sigma := [solve(%)]

[-800, 200, 200]

CharacteristicPolynomial(tau, lambda)

lambda^3+400*lambda^2-280000*lambda+32000000

sigma := sort([solve(%)])

[-800, 200, 200]

NULL

Download mapleac3.mw

@Kitonum But my points were not about how the timing mechanism itself was used, and I don't see how you last comment contradicts any of them.

If you're going to identify numerically computed roots then Roots(...,numeric) isn't much worse of a choice than Analytic(...), and for the example you gave it is slightly faster unless you give one a leg up advantage by adjusting working precision in an ad hoc manner. It doesn't make so much sense to compare against Roots as symbolic solver.

Yes, there are indeed cases in which symbolic Roots(...) "goes away" while attempting to verify symbolically, as you showed. That is the give-and-take between it and symbolic solve. The latter has the problem that it can sometimes miss roots, or (more often) deliver some symbolic results which are not actually roots. We'd all love it is the trig engine within solve were able to compute roots: without need for verification, not miss any, and not return non-roots. It seems difficult.

Verification via shake is a partial solution, but a little tricky to do "right".

First 203 204 205 206 207 208 209 Last Page 205 of 612