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()
Out[61]:
Out[61]:
\[ 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
Out[63]:
\[ x\, D+1 \]

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.