March 26--27, 1998

R. M. Corless

1 (a) Find the fraction-free factorization of

> A := matrix( 3, 3, [2, c, 0, c, 2, 1, 0, 1, 2] );

> fff(A);

The determinant of each F_i is always nonzero, and hence this factorization is always correct. I had intended

to give a more complicated example, but this was ok.

1(b) Consider the eigenvalue problem

> x := vector(3);

> with(linalg):

`Warning, new definition for norm`

`Warning, new definition for trace`

> zero := evalm( A &* x - lambda * x );

> eqs := [seq(zero[i],i=1..3), x[1] - 1];

> with(Groebner):

> gb := gbasis(eqs, tdeg(x[1],x[2],x[3],lambda) );

> wrong := convert( subs( c = 0, gb ), set);

> right := gbasis( subs( c=0, eqs ), tdeg(x[1],x[2],x[3],lambda) );

Question 2.

> interface(echo=2);

`# Gallery(3) matrix`

`> with(linalg):`

`> A := matrix(3,3,[-149,-50,-154,537,180,546,-27,-9,-25]);`

`> E := matrix(3,3,[130,-390,0,43,-129,0,133,-399,0]);`

`> AtE := evalm(A+t*E);`

`> p := charpoly(AtE,mu):`

`> eigs := solve(p,mu):`

> interface(echo=1):

> p := array(1..20):

> pts := map(z->[evalc(Re(z)),evalc(Im(z))],{eigs}):

> for i to 20 do p[i] := plot( subs(t=i/10^7,pts), x=0..4, style=POINT,symbol=CROSS): od:

> with(plots):

> display([seq(p[i],i=1..20)], insequence=true,axes=BOXED);

>

3 (a)

> z := x+I*y;

> f := n -> evalc(Re(z^n));

> with(plots):

> contourplot( f(6), x=-1..1, y=-1..1 , grid=[100,100]);

The others are all similar.

3 (b)

> f := sin( y + sin(x^2*y-1/x) );

>

> contourplot(f, x=-Pi..Pi, y=-Pi..Pi, grid=[100,100],color=black );

This is different from Figure 2.21 on p. 112 of Essential Maple because, although every contour of

> p := y + sin(x^2*y-1/x) = const;

is also a contour of the sine of this, (because sin(const) is also constant), the spacing between the contours

will certainly be different---the contour plotting software chooses constant intervals for the contours, by

default. Also, the contours p and p+2*Pi will be mapped to the same contour in sin(p).

> evalb( sin(P) = sin(P+2*Pi) );

4. Planetary motion.

> restart:

> assume( epsilon > 0, epsilon < 1 );

> r := 1/( 1+epsilon*cos(psi) );

> T := int( r^2, psi=0..theta );

That answer is incorrect, because it contains a spurious discontinuity at theta = Pi/2, whereas the

integrand is continuous everywhere.

> plot( subs(epsilon=1/2, T), theta=0..2*Pi );

> plot3d( T, epsilon=0..1, theta=0..2*Pi, view=[0..1,0..2*Pi,-5..5],
orientation=[25,23]);

> limit(T, theta=Pi, left);

> limit(T, theta=Pi, right);

But the answer differentiates to get the right answer, so we are out only by a piecewise constant.

> simplify( diff( T, theta ) - subs(psi=theta, r^2), trig);

> combine(%);

>