lec26.mws

L-26 MCS 320 Wednesday 26 October 2005

> restart;

The command solve is probably one of the most frequently used commands in Maple.

1. Equations in one variable with a parameter

A strength of computer algebra is that we can deal with parametric problems.

> equ := a^2*(x^2 + x + 1) - a*(2*x - 3) - x^2 -3*x+2;

equ := a^2*(x^2+x+1)-a*(2*x-3)-x^2-3*x+2

This equation is an equation in x (x is the variable), with a as a parameter.

> equ := expand(equ);

equ := a^2*x^2+a^2*x+a^2-2*a*x+3*a-x^2-3*x+2

> sols := solve(equ,x);

sols := -(a-3-(-3*a^2-10*a+17)^(1/2))/(2*(a-1)), -(a-3+(-3*a^2-10*a+17)^(1/2))/(2*(a-1))

The application of the quadratic formula leads to two solutions, as we may expect.  However, this solution is not entirely satisfactory.  For a = 1, the formulas for the solution are no longer valid.

> equ1 := subs(a=1,equ); # we treat the case a = 1 separately

equ1 := -4*x+6

For a = 1 we have the solution

> solve(equ1,x);

3/2

This is the first special case.  In this first special case, we have fewer than two solutions.

As a hint, in search for the second special case, let us organize the equation as an equation in x:

> collect(equ,x);

(a^2-1)*x^2+(a^2-3-2*a)*x+3*a+a^2+2

> subs(a=-1,equ);

0

We discovered -1 as another special value for the parameter a.  In this second special case, the entire equation vanishes.  We are left with 0 = 0.  Any value for x is a solution, i.e.: we have infinitely many solutions.

For this type of problem, in the case of quadratic equations, we can discover the special values for the parameter a automatically:

> solve(equ,a);

-1, (x^2+3*x-2)/(x^2+x+1)

We see the -1 apparing.

2. Difficulties with interpreting the output of solve

In this section we consider more general, nonpolynomial, equations.

the output of solve needs more processing

> eq := x = cos(x);

eq := x = cos(x)

> s := solve(eq,x);

s := RootOf(_Z-cos(_Z))

The RootOf is helpful to declare formal roots of polynomials, but are of no use for general equations.

> evalf(s);

.7390851332

For this type of trigonometic equation, the only outcome which makes sense are numerical approximations.

too few solutions

> eq := sin(x) - 1/2;

eq := sin(x)-1/2

> solve(eq,x);

Pi/6

In this case, Maple ignores the periodicity of the sine function.

> _EnvAllSolutions := true;

_EnvAllSolutions := true

> sols := solve(eq,x);

sols := 1/6*Pi+2/3*Pi*_B1+2*Pi*_Z1

We have a family of solutions, parametrized by the boolean B1 (the underscore before the B indicates that is a variable name chosen by Maple), and by the integer valued parameter Z1.

> family := (b,k) -> subs(_B1=b,_Z1=k,sols);

family := proc (b, k) options operator, arrow; subs(_B1 = b, _Z1 = k, sols) end proc

> family(1,3);

41*Pi/6

too many solutions

> eqn := (sin(x)^3 - 3*sin(x) + 1 = 3);

eqn := sin(x)^3-3*sin(x)+1 = 3

> solve(eqn,x);

arcsin(2)-2*arcsin(2)*_B2+2*Pi*_Z2+Pi*_B2, -1/2*Pi+2*Pi*_Z3, -1/2*Pi+2*Pi*_Z4

But suppose we were only interested in the real solutions.  Then it does not make any sense to ask for the inverse of the sine function at 2.

> evalf(arcsin(2));

1.570796327-1.316957897*I

We should see this equation as a polynomial in sin(x):

> solve(eqn,sin(x));

2, -1, -1

From the output of solve, we ignore the 2 (because we only want real solutions) and keep the -1, which is here counted twice.  This double root -1 was also not obvious from the first output of solve, where we could take different values for Z3 and Z4.

3. Equations in several variables

Here we just illustrate the discovery of Cramer's rule for a system of two linear equations:

> sys := {a[1,1]*x[1]+a[1,2]*x[2] = b[1],
       a[2,1]*x[1]+a[2,2]*x[2] = b[2]};

sys := {a[1, 1]*x[1]+a[1, 2]*x[2] = b[1], a[2, 1]*x[1]+a[2, 2]*x[2] = b[2]}

This is a general 2-by-2 linear system with unknowns {x[1],x[2]}:

> solve(sys,{x[1],x[2]});

{x[1] = (-a[1, 2]*b[2]+a[2, 2]*b[1])/(-a[1, 2]*a[2, 1]+a[2, 2]*a[1, 1]), x[2] = -(-a[1, 1]*b[2]+b[1]*a[2, 1])/(-a[1, 2]*a[2, 1]+a[2, 2]*a[1, 1])}

As before (as in the parametric problem of section 1), the solutions are valid only when the determinant of the matrix is nonzero.  This determinant here appears in the denominator of the solutions.  As before, we distinguish three cases: the general case, the first special case of no solutions, and the second special case of infinitely many equations.

We will return to linear equations when we consider the packages linalg and LinearAlgebra.

4. The Groebner basis method

This is recall of the method we used in connection with Lagrange multipliers.

> p1 := x^2 + x^3 - y^2; p2 := (x-1)^2 + y^2 - 1;

p1 := x^2+x^3-y^2

p2 := (x-1)^2+y^2-1

> sys := [p1,p2];

sys := [x^2+x^3-y^2, (x-1)^2+y^2-1]

> with(Groebner);

[Basis, FGLM, HilbertDimension, HilbertPolynomial, HilbertSeries, InterReduce, IsProper, IsZeroDimensional, LeadingCoefficient, LeadingMonomial, LeadingTerm, MonomialOrder, MultiplicationMatrix, Norma...[Basis, FGLM, HilbertDimension, HilbertPolynomial, HilbertSeries, InterReduce, IsProper, IsZeroDimensional, LeadingCoefficient, LeadingMonomial, LeadingTerm, MonomialOrder, MultiplicationMatrix, Norma...[Basis, FGLM, HilbertDimension, HilbertPolynomial, HilbertSeries, InterReduce, IsProper, IsZeroDimensional, LeadingCoefficient, LeadingMonomial, LeadingTerm, MonomialOrder, MultiplicationMatrix, Norma...[Basis, FGLM, HilbertDimension, HilbertPolynomial, HilbertSeries, InterReduce, IsProper, IsZeroDimensional, LeadingCoefficient, LeadingMonomial, LeadingTerm, MonomialOrder, MultiplicationMatrix, Norma...

> gb := Basis(sys,plex(x,y));

gb := [y^6-12*y^2+12*y^4, 24*x-y^4-18*y^2]

In the pure lexicographical order of the variables in the system, we obtain a triangular system: the first equation is only in y, the second equation in x and y.  We see there are six solutions in complex space.  With solve we can extract the rational solutions first:

> sy := solve(gb[1],y);

sy := 0, 0, -(-6+4*3^(1/2))^(1/2), (-6+4*3^(1/2))^(1/2), -I*(6+4*3^(1/2))^(1/2), (6+4*3^(1/2))^(1/2)*I

We have six values for y.  The 0 occurs twice, so we have a double root at (0,0), look at the second equation.  In addition to the origin, we have two real solutions, let us check

> evalf(sy);

0., 0., -.9634330449, .9634330449, -3.595581070*I, 3.595581070*I

To find the x values, we substitute the values for y into the second equation:

> eq2 := seq(subs(y=sy[k],gb[2]),k=1..nops([sy]));

eq2 := 24*x, 24*x, 24*x-(-6+4*3^(1/2))^2+108-72*3^(1/2), 24*x-(-6+4*3^(1/2))^2+108-72*3^(1/2), 24*x-(6+4*3^(1/2))^2+108+72*3^(1/2), 24*x-(6+4*3^(1/2))^2+108+72*3^(1/2)eq2 := 24*x, 24*x, 24*x-(-6+4*3^(1/2))^2+108-72*3^(1/2), 24*x-(-6+4*3^(1/2))^2+108-72*3^(1/2), 24*x-(6+4*3^(1/2))^2+108+72*3^(1/2), 24*x-(6+4*3^(1/2))^2+108+72*3^(1/2)eq2 := 24*x, 24*x, 24*x-(-6+4*3^(1/2))^2+108-72*3^(1/2), 24*x-(-6+4*3^(1/2))^2+108-72*3^(1/2), 24*x-(6+4*3^(1/2))^2+108+72*3^(1/2), 24*x-(6+4*3^(1/2))^2+108+72*3^(1/2)

> sx := map(eqx->solve(eqx,x),[eq2]);

sx := [0, 0, 3^(1/2)-1, 3^(1/2)-1, -1-3^(1/2), -1-3^(1/2)]

> sols := zip((x,y)->[x,y],sx,[sy]);

sols := [[0, 0], [0, 0], [3^(1/2)-1, -(-6+4*3^(1/2))^(1/2)], [3^(1/2)-1, (-6+4*3^(1/2))^(1/2)], [-1-3^(1/2), -I*(6+4*3^(1/2))^(1/2)], [-1-3^(1/2), (6+4*3^(1/2))^(1/2)*I]]

The solution list sols is a list of lists, with every element in that list the coordinates for x and y.

5. Other solvers

See the lecture note.