Menu
▾
▴
maxima-discuss
|
From: Raymond R. <ray...@gm...> - 2014-02-20 16:59:27
Attachments:
exptest-dump.txt
|
load("simplify_sum")$load(linearalgebra)$
H_Pascal : matrix([0,0,0,0,0],[1,0,0,0,0],[0,2,0,0,0],[0,0,3,0,0],[0,0,0,4,0]);
H_Lag : matrix([0,0,0,0,0],[-1,1,0,0,0],[-1,0,1,0,0],[-1,0,0,1,0],[-1,0,0,0,1]);
P_Pascal : matrixexp(H_Pascal);
P_Lag : matrixexp(H_Lag);
P_Lag_1 : sum((H_Lag^i)/(i!),i,0,4);
|
|
From: Iai M. ax <ia...@ax...> - 2014-03-19 09:24:33
|
Hi, Could you tell me the meaning of this warning message from matrixexp()? "Proviso: assuming 4*(w3^2+w2^2+w1^2) # 0" Can I consider its result is reliable? The commands I entered are below: declare(w1, real); declare(w2, real); declare(w3, real); W:matrix([0, w3, -w2, w1], [-w3, 0, w1, w2], [w2, -w1, 0, w3], [-w1, -w2, -w3, 0]); expW:matrixexp(W); Thanks, Iai |
|
From: Helmut J. <jar...@ig...> - 2014-03-19 11:37:14
|
On 03/19/2014 10:08:12 AM, Iai Masafumi ax wrote: > Hi, > > Could you tell me the meaning of this warning message from > matrixexp()? > > "Proviso: assuming 4*(w3^2+w2^2+w1^2) # 0" > > Can I consider its result is reliable? > > The commands I entered are below: > > declare(w1, real); > declare(w2, real); > declare(w3, real); > W:matrix([0, w3, -w2, w1], [-w3, 0, w1, w2], [w2, -w1, 0, w3], [-w1, > -w2, -w3, 0]); > expW:matrixexp(W); > > factor(determinant(W)); gives (w3^2+w2^2+w1^2)^2, so Maxima assumes that W is a regular matrix. I don't know the internals of matrixexp, but the warning states that the result is OK unless W is singular. Helmut |
|
From: Przemek K. <prz...@ni...> - 2014-03-19 16:05:51
|
On 03/19/2014 05:08 AM, Iai Masafumi ax wrote: > Could you tell me the meaning of this warning message from matrixexp()? > > "Proviso: assuming 4*(w3^2+w2^2+w1^2) # 0" > > Can I consider its result is reliable? > > The commands I entered are below: > > declare(w1, real); > declare(w2, real); > declare(w3, real); > W:matrix([0, w3, -w2, w1], [-w3, 0, w1, w2], [w2, -w1, 0, w3], [-w1, > -w2, -w3, 0]); > expW:matrixexp(W); The expression (2*w3**2+2*w2**2+2*w1**2) appears in the denominator of the solution, so the Proviso warns you that the solution is undefined when 4*(w3^2+w2^2+w1^2) is equal to 0. I agree that the use of # to signify 'unequal to' is quirky---I don't think Maxima uses it as such (I am sure I will be corrected if I am wrong about that). |
|
From: Gunter K. <gu...@pe...> - 2014-03-19 17:30:46
|
Macima actually understand the hash sign (#) as sign for "not equal". Since every programming language uses a different symbol here it seems like there is no standard method in specifying inequality using ASCII characters. On 19. März 2014 17:05:37 MEZ, Przemek Klosowski <prz...@ni...> wrote: >On 03/19/2014 05:08 AM, Iai Masafumi ax wrote: >> Could you tell me the meaning of this warning message from >matrixexp()? >> >> "Proviso: assuming 4*(w3^2+w2^2+w1^2) # 0" >> >> Can I consider its result is reliable? >> >> The commands I entered are below: >> >> declare(w1, real); >> declare(w2, real); >> declare(w3, real); >> W:matrix([0, w3, -w2, w1], [-w3, 0, w1, w2], [w2, -w1, 0, w3], [-w1, >> -w2, -w3, 0]); >> expW:matrixexp(W); >The expression (2*w3**2+2*w2**2+2*w1**2) appears in the denominator of >the solution, so the Proviso warns you that the solution is undefined >when 4*(w3^2+w2^2+w1^2) is equal to 0. I agree that the use of # to >signify 'unequal to' is quirky---I don't think Maxima uses it as such >(I >am sure I will be corrected if I am wrong about that). > > >------------------------------------------------------------------------ > >------------------------------------------------------------------------------ >Learn Graph Databases - Download FREE O'Reilly Book >"Graph Databases" is the definitive new guide to graph databases and >their >applications. Written by three acclaimed leaders in the field, >this first edition is now available. Download your free book today! >http://p.sf.net/sfu/13534_NeoTech > >------------------------------------------------------------------------ > >_______________________________________________ >Maxima-discuss mailing list >Max...@li... >https://lists.sourceforge.net/lists/listinfo/maxima-discuss -- Diese Nachricht wurde von meinem Android-Mobiltelefon mit K-9 Mail gesendet. |
|
From: Barton W. <wi...@un...> - 2014-03-19 17:47:39
|
> Could you tell me the meaning of this warning message from matrixexp()? > "Proviso: assuming 4*(w3^2+w2^2+w1^2) # 0" It means in the course of the calculation there was a division by w3^2+w2^2+w1^2. You might notice that the denominators of exp(W) involve w3^2+w2^2+w1^2. You might like to try this: (%i62) limit(ratsubst(zzz, w3^2+w2^2+w1^2, expW),zzz,0); (%o62) matrix([1,w3,-w2,w1],[-w3,1,w1,w2],[w2,-w1,1,w3],[-w1,-w2,-w3,1]) And I believe that this is correct: When w3^2+w2^2+w1^2, the matrix W is a negative root of unity (W^2 + I = 0). > > Can I consider its result is reliable? > > The commands I entered are below: > > declare(w1, real); > declare(w2, real); > declare(w3, real); > W:matrix([0, w3, -w2, w1], [-w3, 0, w1, w2], [w2, -w1, 0, w3], [-w1, > -w2, -w3, 0]); > expW:matrixexp(W); |
|
From: Iai M. ax <ia...@ax...> - 2014-03-19 22:24:27
|
Thank you guys for replies. It would have been obvious for everybody if the message was Proviso: assuming 4*(w3^2+w2^2+w1^2) not equal to 0 (2014/03/20 2:47), Barton Willis wrote: > >> Could you tell me the meaning of this warning message from matrixexp()? >> "Proviso: assuming 4*(w3^2+w2^2+w1^2) # 0" > > It means in the course of the calculation there was a division by w3^2+w2^2+w1^2. You might > notice that the denominators of exp(W) involve w3^2+w2^2+w1^2. You might like to try this: > > (%i62) limit(ratsubst(zzz, w3^2+w2^2+w1^2, expW),zzz,0); > (%o62) matrix([1,w3,-w2,w1],[-w3,1,w1,w2],[w2,-w1,1,w3],[-w1,-w2,-w3,1]) > > And I believe that this is correct: When w3^2+w2^2+w1^2, the matrix W is a negative root of unity > (W^2 + I = 0). > > >> >> Can I consider its result is reliable? >> >> The commands I entered are below: >> >> declare(w1, real); >> declare(w2, real); >> declare(w3, real); >> W:matrix([0, w3, -w2, w1], [-w3, 0, w1, w2], [w2, -w1, 0, w3], [-w1, >> -w2, -w3, 0]); >> expW:matrixexp(W); |
|
From: Thomas D. D. <to...@wa...> - 2014-03-20 15:37:05
|
I think I have a problem, somewhere. GCL? Maxima version: "5.32.1" Maxima build date: "2014-03-02 18:35:41" Host type: "x86_64-unknown-linux-gnu" Lisp implementation type: "GNU Common Lisp (GCL)" Lisp implementation version: "GCL 2.6.10" kill(all); declare(w1, real), declare(w2, real), declare(w3, real); W:matrix( [ 0, w3, -w2, w1], [-w3, 0, w1, w2], [ w2, -w1, 0, w3], [-w1, -w2, -w3, 0]); expW:matrixexp(W); eq:ratsubst(zzz, w3^2+w2^2+w1^2, expW); limit(eq,zzz,0); and, 12 hours later no result. Core i7-3930k OC 4.2GHz, 16g RAM. When I did ^C, the process was using 1.3g RAM. Tom Dean |
|
From: Richard F. <fa...@be...> - 2014-03-20 16:16:35
|
On 3/20/2014 8:37 AM, Thomas D. Dean wrote: > I think I have a problem, somewhere. GCL? > > Maxima version: "5.32.1" > Maxima build date: "2014-03-02 18:35:41" > Host type: "x86_64-unknown-linux-gnu" > Lisp implementation type: "GNU Common Lisp (GCL)" > Lisp implementation version: "GCL 2.6.10" > > kill(all); > declare(w1, real), declare(w2, real), declare(w3, real); > W:matrix( > [ 0, w3, -w2, w1], > [-w3, 0, w1, w2], > [ w2, -w1, 0, w3], > [-w1, -w2, -w3, 0]); > > expW:matrixexp(W); > eq:ratsubst(zzz, w3^2+w2^2+w1^2, expW); > limit(eq,zzz,0); > > and, 12 hours later no result. > > Core i7-3930k OC 4.2GHz, 16g RAM. > > When I did ^C, the process was using 1.3g RAM. > > Tom Dean > > If we assume that as zzz -> 0, so does sqt(-zzz), then perhaps you could just do ratsubst(yyy, sqrt(-zzz),eq). then limit(%,yyy,0) returns the 4x4 identity matrix in a fraction of a second. |
|
From: Richard F. <fa...@be...> - 2014-03-20 17:56:42
|
On 3/20/2014 9:16 AM, Richard Fateman wrote: > > If we assume that as zzz -> 0, so does sqt(-zzz), then perhaps you could > just do ratsubst(yyy, sqrt(-zzz),eq). > > then limit(%,yyy,0) returns the 4x4 identity matrix in a fraction of a > second. oops, I see that the ratsubst(yyy, sqrt(-zzz) , eq) left some instances of zzz around, so the limits may not be right. you can play around some more with something like *for i:1 thru 4 do for j:1 thru 4 do print([i,j,limit(m[i,j],yyy,0)]);* |
|
From: Robert D. <rob...@gm...> - 2014-03-20 18:25:38
|
On 2014-03-20, Thomas D. Dean <to...@wa...> wrote:
> declare(w1, real), declare(w2, real), declare(w3, real);
> W:matrix(
> [ 0, w3, -w2, w1],
> [-w3, 0, w1, w2],
> [ w2, -w1, 0, w3],
> [-w1, -w2, -w3, 0]);
>
> expW:matrixexp(W);
> eq:ratsubst(zzz, w3^2+w2^2+w1^2, expW);
> limit(eq,zzz,0);
Looks like the problem is in limit. Here's a simpler version of it:
(%i1) foo : sqrt(-z)*%e^-sqrt(-z)/z-sqrt(-z)*%e^sqrt(-z)/z $
(%i2) limit (foo, z, 0);
<waits indefinitely>
Tracing the internal function LIMIT, I see that LIMIT is called on
enormous expressions which contain large integers. I'm not sure where
those come from. Factoring something? Rationalizing bigfloat
coefficients?
Incidentally gruntz(foo, z, 0, plus) => 2.
At this point all I can suggest is that someone file a bug report.
best,
Robert Dodier
|
|
From: Robert D. <rob...@gm...> - 2014-02-20 19:47:19
|
On 2014-02-20, Raymond Rogers <ray...@gm...> wrote:
> load("simplify_sum")$load(linearalgebra)$
> H_Pascal : matrix([0,0,0,0,0],[1,0,0,0,0],[0,2,0,0,0],[0,0,3,0,0],[0,0,0,4,0]);
> H_Lag : matrix([0,0,0,0,0],[-1,1,0,0,0],[-1,0,1,0,0],[-1,0,0,1,0],[-1,0,0,0,1]);
> P_Pascal : matrixexp(H_Pascal);
> P_Lag : matrixexp(H_Lag);
> P_Lag_1 : sum((H_Lag^i)/(i!),i,0,4);
Raymond, don't you mean to write H_Lag^^i (matrix power) in the
summation? I find that when I write it that way, I get a result
very close to what matrixexp returns.
best,
Robert Dodier
|
|
From: Mark <zei...@ya...> - 2014-02-20 22:05:10
|
On 02/20/2014 05:59 PM, Raymond Rogers wrote: > Because of the way matexp is implemented the following produces the > wrong result: > > The following has an incorrect result > Two files are attachedon: e source and one window dump generated from > the source. > H_Lag : > matrix([0,0,0,0,0],[-1,1,0,0,0],[-1,0,1,0,0],[-1,0,0,1,0],[-1,0,0,0,1]); > P_Lag : matrixexp(H_Lag); > > [ 1 0 0 0 0 ] > [ ] > [ 1 - %e %e 0 0 0 ] > [ ] > [ 1 - %e 0 %e 0 0 ] > [ ] > [ 1 - %e 0 0 %e 0 ] > [ ] > [ 1 - %e 0 0 0 %e ] > > Which is certainly wrong Are you certain? Your matrix H_Lag satisfies H_Lag^^n = H_Lag for any integer n>0. So, matrixexp(H_Lag) = identfor(H) + H_Lag*(%e-1) which is your P_Lag -M |
|
From: Barton W. <wi...@un...> - 2014-02-21 09:29:34
|
(Initially I didn't reply to all.) (%i51) M : matrix([0,0,0,0,0],[-1,1,0,0,0],[-1,0,1,0,0],[-1,0,0,1,0],[-1,0,0,0,1])$ The matrix M is a projection; let's check: (%i52) every(lambda([s], is(s=0)), M.M- M); (%o52) true So exp(M) = I/0! + M/1! + M^2/2! + ... = I + (1/1! + 1/2! + ...) * M = I (%e - 1) M. Let's compare this to matrixexp(M) (%i53) identfor(M) + (%e-1) * M - matrixexp(M); (%o53) matrix([0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0]) This looks OK to me. What am I missing? --Barton (author of matrixexp) |
|
From: Raymond R. <ray...@gm...> - 2014-02-21 02:24:43
|
On 02/20/2014 05:05 PM, Mark wrote: > On 02/20/2014 05:59 PM, Raymond Rogers wrote: >> Because of the way matexp is implemented the following produces the >> wrong result: >> >> The following has an incorrect result >> Two files are attachedon: e source and one window dump generated from >> the source. >> H_Lag : >> matrix([0,0,0,0,0],[-1,1,0,0,0],[-1,0,1,0,0],[-1,0,0,1,0],[-1,0,0,0,1]); >> P_Lag : matrixexp(H_Lag); >> >> [ 1 0 0 0 0 ] >> [ ] >> [ 1 - %e %e 0 0 0 ] >> [ ] >> [ 1 - %e 0 %e 0 0 ] >> [ ] >> [ 1 - %e 0 0 %e 0 ] >> [ ] >> [ 1 - %e 0 0 0 %e ] >> >> Which is certainly wrong > Are you certain? i am absolutely certain I made a mistake and apologize for wasting peoples time thinking about it! Embarrassed Ray > > Your matrix H_Lag satisfies H_Lag^^n = H_Lag for any integer n>0. So, > matrixexp(H_Lag) = identfor(H) + H_Lag*(%e-1) > which is your P_Lag > > -M -- Act IV, Sc. IV What is a man, If his chief good and profit of his time Be to sleep and feed. Be a beast, no more |
|
From: Barton W. <wi...@un...> - 2014-02-22 15:12:17
|
Incidentally, matrixexp works by decomposing matrices as a linear combination of commuting projections plus a nilpotent. The projections are the (matrix valued) residues at each simple pole of the resolvent matrix. Demonstrating that these residues are commuting projections is a comfort math exercise involving facts about the partial fraction decomposition. This method is based on the linear algebra review in the first chapter of Perturbation Theory for Linear Operators (Tosio Kato). The treatment in Kato uses the Cauchy integral theorem applied to the resolvent instead of the partial fraction decomposition. I do not recommend using matrixexp with floating point entries. There are of course, many other methods for computing (or approximating) the matrix exponential. In part, I wrote matrixexp because I enjoy the treatment in the (now) Classic text by Kato. --Barton |
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X