FriCAS-Solve.input

-- ---
-- jupyter:
--   jupytext:
--     formats: ipynb,input:light
--     text_representation:
--       extension: .input
--       format_name: light
--       format_version: '1.5'
--       jupytext_version: 1.10.0
--   kernelspec:
--     display_name: FriCAS
--     language: spad
--     name: jfricas
-- ---

-- This notebook is licenced under
-- [CC BY-SA 4.0](http://creativecommons.org/licenses/by-sa/4.0/).
--
-- # FriCAS Tutorial (Solve)
--
-- ## Ralf Hemmecke <[ralf@hemmecke.org](mailto:ralf@hemmecke.org)>
--
-- Sources at [Github](http://github.com/fricas/fricas-notebooks/).

)read init.input )quiet

-- # Solving Systems of Polynomial Equations

p1 := x^2 + y^2 + z^2 - 9;
p2 := (x-1)^2 + (y -2)^2 - z;
p3 := y - 1;
lp := [p1, p2, p3]

-- The form of the solution depends on the order of the variables.

solve(lp,[x,y,z])

solve([p1,p2,p3],[z,y,x])

-- A Gröbner basis computation is done with respect to a
-- lexicographical order of all variables that appear in the expression.

gb := groebner(lp)

[factor p for p in gb]

-- 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.

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

-- 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.

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

-- FriCAS also implements a Gröbner factorizer.
-- That is an algorithm which factorizes polynomials during a
-- Gröbner basis computation.

groebnerFactorize(lpl)

-- ## Operation on polynomial ideals

Q := Fraction Integer;
MP := GeneralDistributedMultivariatePolynomial(vl, Q, ET);
V := OrderedVariableList vl;
PId := PolynomialIdeal(Q, ET, V, MP);

mi: List(MP) := [x^2 + y^2 + z^2 - 9, (x-1)^2 + (y -2)^2 - z]
I := mi :: PId

dimension I

mj: List(MP) := [y^2-1, x*y -z]
J := mj :: PId

dimension J

-- The sum of two ideals $I$ and $J$ is the smallest ideal of the ring,
-- that contains $I$ and $J$.

K := I+J

dimension K

-- The [product of two polynomial ideals](
-- http://en.wikipedia.org/wiki/Ideal_%28ring_theory%29#Ideal_operations)
-- can be computed as well.

K := I * J

L := intersect(I, J)

(K = L)@Boolean

-- The leading ideal of a polynomial ideal $K$ is the ideal that is
-- generated by all leading terms of elements of $K$.

leadingIdeal K

leadingIdeal L

-- 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$.

f: MP := y-1
saturate(J,f)

-- Of course, if we saturate a polynomial ideal $J$ by an element
-- of $J$, we get the whole ring.

saturate(J, mj.1)

-- ## Solving Diophantine Equations

m := matrix [[1,2],[-1,3]]
v: Vector(Integer) := vector [8, 7]
diophantineSystem(m, v)

-- ## Solving Ordinary Differential Equations

-- We first must tell FriCAS that `y` should be considered as
-- an unknown function.

y := operator 'y

deq := D(y(x), x, 2) + y(x)

-- There are, in fact, infinitely many solutions to this equation,
-- so FriCAS returns a particular solution together with basis
-- for the corresponding homogeneous equation.

solve(deq, y, x)

-- 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$.

solve(deq, y, x=0, [1,2])

-- The coefficients of the differential equation can also be
-- rational functions.
--
-- To make checking later easier, we define a function $d$.

d(f) == x/2*D(f,x,2) - D(f,x) + (x^2+1)/x * f = x^2

deq := d(y(x))

sol := solve(deq, y, x)

-- 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.

E := Expression(Integer);
s := sol :: Record(particular: E, basis: List(E));
g(a,b) == s.particular + a*s.basis.1 + b*s.basis.2

g(a,b)

d g(a,b)

-- The above function $d$ looks like a differential operator,
-- but its type tells something else.

d

-- Working with true differential operators is possible.

Q := Fraction Integer
P := UnivariatePolynomial('x, Q)

-- The polynomial ring $P$ has a derivation `D: % -> %`.

L := LinearOrdinaryDifferentialOperator(P, D)
Dx: L := D()

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

L has Ring

-- 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.

Dx*x

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

Dx(x)

-- 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.