>From mailnull Wed Jan 2 19:05:12 2013
From: Rupert Swarbrick <rswarbrick at gmail.com>
Date: Wed, 2 Jan 2013 18:53:12 +0000
Content-Type: multipart/mixed;
boundary="===============0252527135626607063=="
--===============0252527135626607063==
Content-Type: multipart/signed; boundary="=-=-="; micalg=pgp-sha1;
protocol="application/pgp-signature"
--=-=-=
Content-Type: text/plain
Content-Transfer-Encoding: quoted-printable
Evan Cooch <evan.cooch at gmail.com> writes:
> Greetings --
>
> Consider the following matrix:
>
> a : matrix([-0.7-0.3*h1,0.6-0.6*h1],[0.5-0.5*h2,-0.35-0.65*h2]);
>
> I'm trying to solve for h1 and h2, such that dominant eigenvalue of a is =
0.
>
> I know how to do this in a couple of other CAS systems, but am stumped
> with Maxima. I know the solution in thiscase is
>
> h2=3D(11-81*h1)/(151-21*h1)
>
> Any pointers would be much appreciated (to aid and abet my attempt to
> port a bunch of teaching material to Maxima).
>
> Thanks!
(Apologies for multiple posts)
This is somewhat more involved than your solution suggests. You can
also do everything using just the core Maxima, without to_poly_solve.
(%i1) display2d:false$
(%i2) a : rat(matrix([-0.7-0.3*h1,0.6-0.6*h1],[0.5-0.5*h2,-0.35-0.65*h2]))$
(%i3) p:charpoly(a,t);
(%o3) (200*t^2+(130*h2+60*h1+210)*t+(-21*h1+151)*h2+81*h1-11)/200
(%i4) r:solve(subst(t=0,p),h2);
(%o4) [h2 = (81*h1-11)/(21*h1-151)]
This gives a relationship between h2 & h1 when det(a)=0, i.e. a has 0 as an eigenvalue.
If 0 is the dominant eigenvalue, then you must look at when the 2nd ev is negative:
(%i5) z:rat(subst(r,p));
(%o5) ((210*h1-1510)*t^2+(63*h1^2+294*h1-1657)*t)/(210*h1-1510)
(%i6) e:solve(z,t);
(%o6) [t = -(63*h1^2+294*h1-1657)/(210*h1-1510),t = 0]
You see the non-zero eigenvalue is a rational function of h1; you need
to see when it is negative.
(%i7) fourier_elim([rhs(first(e)) < 0],[h1]);
(%o7) [151/21 < h1,63*h1^2+294*h1-1657 > 0]
or [h1 < 151/21,-(63*h1^2+294*h1-1657) > 0]
By my calculation, this is a union of two open intervals,
(151/21,infty)
(r1,r2)
where r1 < 0 < r2 are the roots of the quadratic.
Leo
An almost complete answer:
Firstly, let's not use the icky decimals any longer than necessary:
(%i1) a : rat(matrix([-0.7-0.3*h1,0.6-0.6*h1],[0.5-0.5*h2,-0.35-0.65*h2]));
<snip: warnings from rat>
[ 3 h1 + 7 3 h1 - 3 ]
[ - -------- - -------- ]
[ 10 5 ]
(%o1)/R/ [ ]
[ h2 - 1 13 h2 + 7 ]
[ - ------ - --------- ]
[ 2 20 ]
Right, now use eivals to grab the eigenvalues:
(%i2) load("eigen")$
(%i3) evals: first(eivals(a));
2 2
(%o3) [- (sqrt(169 h2 + (324 h1 - 662) h2 + 36 h1 - 396 h1 + 529) + 13 h2
+ 6 h1 + 21)/40,=20
2 2
sqrt(169 h2 + (324 h1 - 662) h2 + 36 h1 - 396 h1 + 529) - 13 h2 - 6 h1 - =
21
=2D------------------------------------------------------------------------=
----]
40
The solve() function doesn't seem to have much luck, but to_poly_solve
does much better:
(%i4) load("to_poly_solve")$
*** output flushed ***
(%i5) to_poly_solve(first(evals)=3D0, h2);
2
%pi 126 h1 + 588 h1 - 3314
(%o5) %union(%if((- --- < parg(- -----------------------))
2 21 h1 - 151
2
126 h1 + 588 h1 - 3314 %pi 81 h1 - 11
%and (parg(- -----------------------) <=3D ---), [h2 =3D -----------], %u=
nion()))
21 h1 - 151 2 21 h1 - 151
(%i6) to_poly_solve(second(evals)=3D0, h2);
2
%pi 126 h1 + 588 h1 - 3314
(%o6) %union(%if((- --- < parg(-----------------------))
2 21 h1 - 151
2
126 h1 + 588 h1 - 3314 %pi 81 h1 - 11
%and (parg(-----------------------) <=3D ---), [h2 =3D -----------], %u=
nion()))
21 h1 - 151 2 21 h1 - 151
Frustratingly, I can't see how to tell Maxima that h1 and h2 are real,
so all the parg values are zero. Also, the %union with only one element
is a bit rubbish. Anyway, I can do one more thing:
(%i7) eq: first(second(first(%o6)));
81 h1 - 11
(%o7) h2 =3D -----------
21 h1 - 151
Then:
(%i21) radcan(second(evals)), eq;
(%o21) 0
Looking at the picture given by
plot2d(first(evals), [h1, -10, 5]), eq;
it looks like something more complicated is going on with the other
eigenvalue. Maybe you know a bit more about the situation - I'd be
interested to know what the graph is saying.
Rupert
--=-=-=
Content-Type: application/pgp-signature
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.12 (GNU/Linux)
iJwEAQECAAYFAlDkghgACgkQRtd/pJbYVoanwAP/Tf356BZKOeJH8XIroEjtkdIu
2W6V0g/Rf2LAKaCN5IEditQQ133/W94Ljt7KtW41cWS/EcJcFX5O+SU/nosLvEDI
KQAqYIX55MMmbb4O9pZ9WZNOyQMDwpMP9wlyrg2vtY2I6I3lnRHLptl3cV0ZT1JI
pp7jt4CKYEeLfnVOLg8=
=0Mm4
-----END PGP SIGNATURE-----
--=-=-=--
--===============0252527135626607063==
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline
_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
--===============0252527135626607063==--