This notebook is licenced under CC BY-SA 4.0.

# FriCAS Tutorial (Solve)¶

## Ralf Hemmecke <ralf@hemmecke.org>¶

Sources at Github.

In [1]:
)set message type off
)set output algebra off
setFormat!(FormatMathJax)$JFriCASSupport )set message type on In [2]: )version Value = "FriCAS a9422d32eb6f6b03d98ac6adfc9bf05ec0f5e8bd compiled at Sa 03 Apr 2021 00:29:13 CEST" # Solving Systems of Polynomial Equations¶ In [5]: p1 := x^2 + y^2 + z^2 - 9; p2 := (x-1)^2 + (y -2)^2 - z; p3 := y - 1; lp := [p1, p2, p3] Out[5]: Out[5]: Out[5]: Out[5]: $\left[{z}^{2}+{y}^{2}+{x}^{2}-9, -z+{y}^{2}-4\, y+{x}^{2}-2\, x+5, y-1\right]$ The form of the solution depends on the order of the variables. In [6]: solve(lp,[x,y,z]) Out[6]: $\left[\left[x=2, y=1, z=2\right], \left[x=\frac{-{z}^{2}-z+10}{2}, y=1, {z}^{3}+4\, {z}^{2}-7\, z-34=0\right]\right]$ In [7]: solve([p1,p2,p3],[z,y,x]) Out[7]: $\left[\left[z=2, y=1, x=2\right], \left[z={x}^{2}-2\, x+2, y=1, {x}^{3}-2\, {x}^{2}+5\, x+2=0\right]\right]$ A Gröbner basis computation is done with respect to a lexicographical order of all variables that appear in the expression. In [8]: gb := groebner(lp) Out[8]: $\left[z-{x}^{2}+2\, x-2, y-1, {x}^{4}-4\, {x}^{3}+9\, {x}^{2}-8\, x-4\right]$ In [9]: [factor p for p in gb] Out[9]: $\left[z-{x}^{2}+2\, x-2, y-1, \left(x-2\right)\, \left({x}^{3}-2\, {x}^{2}+5\, x+2\right)\right]$ Depending on the type of the polynomials the Groebner package is able to compute with respect to a lexicographical termorder as specified by the exponent domain EL. In [15]: vl := [x,y,z]; EL := DirectProduct(#vl, NonNegativeInteger); PL := GeneralDistributedMultivariatePolynomial(vl, Integer, EL); GL := GroebnerPackage(Integer, EL, PL); lpl: List(PL) := [x^2 + y^2 + z^2 - 9, (x-1)^2 + (y -2)^2 - z, y-1]; groebner(lpl)$GL
Out[15]:
Out[15]:
Out[15]:
Out[15]:
Out[15]:
Out[15]:
$\left[2\, x+{z}^{2}+z-10, y-1, {z}^{4}+2\, {z}^{3}-15\, {z}^{2}-20\, z+68\right]$

A Gröbner basis computation with respect to a degree reverse lexicographical order just mean to provide an exponent domain that implements the right term order.

In [20]:
ET := HomogeneousDirectProduct(#vl, NonNegativeInteger);
PT := GeneralDistributedMultivariatePolynomial(vl, Integer, ET);
GT := GroebnerPackage(Integer, ET, PT);
lpt: List(PT) := [x^2 + y^2 + z^2 - 9, (x-1)^2 + (y -2)^2 - z, y-1];
groebner(lpt)$GT Out[20]: Out[20]: Out[20]: Out[20]: Out[20]: $\left[{x}^{2}-2\, x-z+2, {z}^{2}+2\, x+z-10, y-1\right]$ FriCAS also implements a Gröbner factorizer. That is an algorithm which factorizes polynomials during a Gröbner basis computation. In [21]: groebnerFactorize(lpl) Out[21]: $\left[\left[2\, x+{z}^{2}+z-10, y-1, {z}^{3}+4\, {z}^{2}-7\, z-34\right], \left[x-2, y-1, z-2\right]\right]$ ## Operation on polynomial ideals¶ In [25]: Q := Fraction Integer; MP := GeneralDistributedMultivariatePolynomial(vl, Q, ET); V := OrderedVariableList vl; PId := PolynomialIdeal(Q, ET, V, MP); Out[25]: Out[25]: Out[25]: Out[25]: In [27]: mi: List(MP) := [x^2 + y^2 + z^2 - 9, (x-1)^2 + (y -2)^2 - z] I := mi :: PId Out[27]: $\left[{x}^{2}+{y}^{2}+{z}^{2}-9, {x}^{2}+{y}^{2}-2\, x-4\, y-z+5\right]$ Out[27]: $\left[{x}^{2}+{y}^{2}+{z}^{2}-9, {x}^{2}+{y}^{2}-2\, x-4\, y-z+5\right]$ In [28]: dimension I Out[28]: $1$ In [30]: mj: List(MP) := [y^2-1, x*y -z] J := mj :: PId Out[30]: $\left[{y}^{2}-1, x\, y-z\right]$ Out[30]: $\left[{y}^{2}-1, x\, y-z\right]$ In [31]: dimension J Out[31]: $1$ The sum of two ideals$I$and$J$is the smallest ideal of the ring, that contains$I$and$J$. In [32]: K := I+J Out[32]: $\left[x-2, y-1, z-2\right]$ In [33]: dimension K Out[33]: $0$ The product of two polynomial ideals can be computed as well. In [34]: K := I * J Out[34]: $\left[{x}^{4}-{y}^{4}-{z}^{4}+4\, x\, {y}^{2}+8\, {y}^{3}+2\, {y}^{2}\, z-8\, {x}^{2}-18\, {y}^{2}+8\, {z}^{2}-4\, x-8\, y-2\, z+19, {x}^{3}\, y+x\, {y}^{3}-2\, {x}^{2}\, y-4\, x\, {y}^{2}-{x}^{2}\, z-x\, y\, z-{y}^{2}\, z+5\, x\, y+2\, x\, z+4\, y\, z+{z}^{2}-5\, z, {x}^{2}\, {y}^{2}+{y}^{4}-2\, x\, {y}^{2}-4\, {y}^{3}-{y}^{2}\, z-{x}^{2}+4\, {y}^{2}+2\, x+4\, y+z-5, {x}^{3}\, z+x\, {y}^{2}\, z-x\, {z}^{3}+{x}^{3}-14\, {x}^{2}\, y+x\, {y}^{2}-14\, {y}^{3}-2\, {x}^{2}\, z-8\, x\, y\, z+2\, {y}^{2}\, z-x\, {z}^{2}-6\, y\, {z}^{2}+2\, {z}^{3}+4\, {x}^{2}+16\, x\, y+36\, {y}^{2}+19\, x\, z+8\, y\, z+4\, {z}^{2}-9\, x+14\, y-18\, z-36, {x}^{2}\, y\, z+{y}^{3}\, z-{x}^{3}-x\, {y}^{2}-2\, x\, y\, z-4\, {y}^{2}\, z-y\, {z}^{2}+2\, {x}^{2}+4\, x\, y+x\, z+5\, y\, z-5\, x, {x}^{2}\, {z}^{2}-{z}^{4}+2\, {x}^{3}+4\, {x}^{2}\, y+{x}^{2}\, z-2\, x\, {z}^{2}-4\, y\, {z}^{2}-{z}^{3}-14\, {x}^{2}+14\, {z}^{2}, x\, y\, {z}^{2}+2\, {x}^{2}\, y+4\, x\, {y}^{2}+x\, y\, z-{z}^{3}-14\, x\, y-2\, x\, z-4\, y\, z-{z}^{2}+14\, z, {y}^{2}\, {z}^{2}+2\, x\, {y}^{2}+4\, {y}^{3}+{y}^{2}\, z-14\, {y}^{2}-{z}^{2}-2\, x-4\, y-z+14, y\, {z}^{3}+2\, x\, y\, z+4\, {y}^{2}\, z-x\, {z}^{2}+y\, {z}^{2}-2\, {x}^{2}-4\, x\, y-x\, z-14\, y\, z+14\, x\right]$ In [35]: L := intersect(I, J) Out[35]: $\left[{x}^{2}\, {y}^{2}+{y}^{4}-2\, x\, {y}^{2}-4\, {y}^{3}-{y}^{2}\, z-{x}^{2}+4\, {y}^{2}+2\, x+4\, y+z-5, {x}^{2}\, y\, z+{y}^{3}\, z+5\, {x}^{2}\, y+5\, {y}^{3}+{x}^{2}\, z-2\, x\, y\, z-3\, {y}^{2}\, z-x\, {z}^{2}-{z}^{3}+3\, {x}^{2}-12\, x\, y-11\, {y}^{2}-5\, x\, z-7\, y\, z-{z}^{2}+6\, x-5\, y+15\, z+11, {x}^{2}\, {z}^{2}-{z}^{4}-6\, {x}^{2}\, y-2\, x\, {y}^{2}-10\, {y}^{3}-{x}^{2}\, z-2\, {y}^{2}\, z-6\, y\, {z}^{2}+{z}^{3}-16\, {x}^{2}+32\, x\, y+22\, {y}^{2}+12\, x\, z+24\, y\, z+16\, {z}^{2}-22\, x+10\, y-30\, z-22, x\, y\, {z}^{2}+2\, {x}^{2}\, y+4\, x\, {y}^{2}+x\, y\, z-{z}^{3}-14\, x\, y-2\, x\, z-4\, y\, z-{z}^{2}+14\, z, {y}^{2}\, {z}^{2}+2\, x\, {y}^{2}+4\, {y}^{3}+{y}^{2}\, z-14\, {y}^{2}-{z}^{2}-2\, x-4\, y-z+14, y\, {z}^{3}+2\, x\, y\, z+4\, {y}^{2}\, z-x\, {z}^{2}+y\, {z}^{2}-2\, {x}^{2}-4\, x\, y-x\, z-14\, y\, z+14\, x, {x}^{3}+5\, {x}^{2}\, y+x\, {y}^{2}+5\, {y}^{3}+{x}^{2}\, z+{y}^{2}\, z-x\, {z}^{2}+y\, {z}^{2}-{z}^{3}+{x}^{2}-16\, x\, y-11\, {y}^{2}-6\, x\, z-12\, y\, z-{z}^{2}+11\, x-5\, y+15\, z+11\right]$ In [36]: (K = L)@Boolean Out[36]: $\texttt{false}$ The leading ideal of a polynomial ideal$K$is the ideal that is generated by all leading terms of elements of$K$. In [37]: leadingIdeal K Out[37]: $\left[{x}^{4}, {x}^{3}\, y, {x}^{2}\, {y}^{2}, {x}^{3}\, z, {x}^{2}\, y\, z, {x}^{2}\, {z}^{2}, x\, y\, {z}^{2}, {y}^{2}\, {z}^{2}, y\, {z}^{3}\right]$ In [38]: leadingIdeal L Out[38]: $\left[{x}^{2}\, {y}^{2}, {x}^{2}\, y\, z, {x}^{2}\, {z}^{2}, x\, y\, {z}^{2}, {y}^{2}\, {z}^{2}, y\, {z}^{3}, {x}^{3}\right]$ The ideal saturate(J, f) contains every polynomial$p$of the underlying polynomial ring MP for which there exists a natural number$n$such that$(f^n)\cdot p \in J$. In [40]: f: MP := y-1 saturate(J,f) Out[40]: $y-1$ Out[40]: $\left[x+z, y+1\right]$ Of course, if we saturate a polynomial ideal$J$by an element of$J$, we get the whole ring. In [41]: saturate(J, mj.1) Out[41]: $\left[1\right]$ ## Solving Diophantine Equations¶ In [44]: m := matrix [[1,2],[-1,3]] v: Vector(Integer) := vector [8, 7] diophantineSystem(m, v) Out[44]: $\begin{bmatrix}1&2\\-1&3\end{bmatrix}$ Out[44]: $\left[8, 7\right]$ Out[44]: $\left[particular=\left[2, 3\right], basis=\left[\right]\right]$ Record(particular: Union(Vector(Integer),"failed"),basis: List(Vector(Integer))) ## Solving Ordinary Differential Equations¶ We first must tell FriCAS that y should be considered as an unknown function. In [45]: y := operator 'y Out[45]: $y$ In [46]: deq := D(y(x), x, 2) + y(x) Out[46]: $y''\left(x\right)+y\left(x\right)$ There are, in fact, infinitely many solutions to this equation, so FriCAS returns a particular solution together with basis for the corresponding homogeneous equation. In [47]: solve(deq, y, x) Out[47]: $\left[particular=0, basis=\left[\cos\left(x\right), \sin\left(x\right)\right]\right]$ Union(Record(particular: Expression(Integer),basis: List(Expression(Integer))),...) Giving initial conditions allows to select a certain solution. For example, we want a solution$y(x)$so that$y(0)=1$and$y'(0)=2$. In [48]: solve(deq, y, x=0, [1,2]) Out[48]: $2\, \sin\left(x\right)+\cos\left(x\right)$ The coefficients of the differential equation can also be rational functions. To make checking later easier, we define a function$d$. In [49]: d(f) == x/2*D(f,x,2) - D(f,x) + (x^2+1)/x * f = x^2 Out[49]: In [50]: deq := d(y(x)) Compiling function d with type Expression(Integer) -> Equation(Expression( Integer)) Out[50]: $\frac{{x}^{2}\, y''\left(x\right)-2\, x\, y'\left(x\right)+\left(2\, {x}^{2}+2\right)\, y\left(x\right)}{2\, x}={x}^{2}$ In [51]: sol := solve(deq, y, x) Out[51]: $\left[particular=x, basis=\left[x\, {e}^{x\, \sqrt{-2}}, x\, {e}^{-2\, x\, \sqrt{-2}}\, {e}^{x\, \sqrt{-2}}\right]\right]$ Union(Record(particular: Expression(Integer),basis: List(Expression(Integer))),...) Let's check the result. Since the result is of type Union, we first have to get rid of it so that the definition of the generic solution$g(a,b)$type checks. In [54]: E := Expression(Integer); s := sol :: Record(particular: E, basis: List(E)); g(a,b) == s.particular + a*s.basis.1 + b*s.basis.2 Out[54]: Out[54]: Out[54]: In [55]: g(a,b) Compiling function g with type (Variable(a), Variable(b)) -> Expression( Integer) Out[55]: $\left(b\, x\, {e}^{-2\, x\, \sqrt{-2}}+a\, x\right)\, {e}^{x\, \sqrt{-2}}+x$ In [56]: d g(a,b) Out[56]: ${x}^{2}={x}^{2}$ The above function$d$looks like a differential operator, but its type tells something else. In [57]: d Out[57]: $d\ f\ ==\ \frac{x}{2}\, D\left(f, x, 2\right)-D\left(f, x\right)+\frac{{x}^{2}+1}{x}\, f={x}^{2}$ Working with true differential operators is possible. In [59]: Q := Fraction Integer P := UnivariatePolynomial('x, Q) Out[59]: Out[59]: The polynomial ring$P\$ has a derivation D: % -> %.

In [61]:
L := LinearOrdinaryDifferentialOperator(P, D)
Dx: L := D()

Of course, this is a ring, but a non-commutative one.

In [62]:
L has Ring
Out[62]:
$\texttt{true}$

The multiplication sign is again overloaded. It's meaning depends on its operands. Internally, linear diffferential operators are stored in a canonical form with coefficients appearing on the left of the D operator.

In [63]:
Dx*x

Of course, applying an operator is not the same as multiplying two operators. Look at the result type.

In [64]:
Dx(x)
Out[64]:
$1$

It should be clear that one can construct linear differential operators not only over a polynomial ring, but also over rational functions. Even operators with matrix coefficient or power series coefficients work.