DMelt:Symbolic/3 JScl
Using the JScl engine
The symbolic capabilities of the "jscl" engine are currently:
- polynomial system solving
- vectors and matrices
- factorization
- derivatives
- integrals (rational functions)
- boolean algebra
- simplification
- geometric algebra
- java code generation
- graphical rendering of math expressions
Symbolic calculations using "jscl" engine is initialized as
from jhplot.math import *
j=Symbolic("jscl") # this sets the engine to "jscl" (R.Jolly etc.)
It allows vector and matrix expressions processing. Vectors are noted as {a,b,c} and matrices as {{a,b},{c,d}}. The usual function is one of:
sin cos tan cot asin acos atan acot log exp sqrt cubic sinh cosh tanh coth asinh acosh atanh acoth abs sgn conjugate eq le ge ne lt gt ap
A reference
| Constants | pi=3.141592653589793 |
| Variables | a, ai a[i1]...[in] a', a'...' |
| Vectors and matrices | {a1,...,an} {{a11,...,a1p},...,{an1,...,anp}} |
| Expressions | a+b, a-b, a*b, a/b, ab, a! div(a,b), mod(a,b) a*(b+c) |
| Trigonometric and hyperbolic functions | sin(x), cos(x), tan(x), cot(x) asin(x), acos(x), atan(x), acot(x) sinh(x), cosh(x), tanh(x), coth(x) asinh(x), acosh(x), atanh(x), acoth(x) |
| Logarithm and exponential | log(x), exp(x) |
| Absolute value, sign | abs(x), sgn(x) |
| Square root, cubic root | sqrt(x), cubic(x) |
| Conjugate | conjugate(x) |
| Implicit roots | root[i](a0,...,an) |
| Implicit functions | f(x), f'(x), f'...'(x) f(x1,...,xm), f{n1,...,nm}(x1,...,xm) f[i1]...[in](x) |
| Derivation | d(f(x),x,value_opt,order_opt) |
| Integration | integral(f(x),x,a,b), integral(f(x),x) |
Vector algebra
grad(f(x1,...,xn),{x1,...,xn}) diverg({f1(x1,...,xn),...,fn(x1,...,xn)},{x1,...,xn}) curl({fx(x,y,z),fy(x,y,z),fz(x,y,z)},{x,y,z})
laplacian(f(x1,...,xn),{x1,...,xn})
dalembertian(f(t,x,y,z),{t,x,y,z})
jacobian({f1(x1,...,xn),...,fm(x1,...,xn)},{x1,...,xn})
Scalar and vector product
{a1,...,an}*{b1,...,bn} vector({ax,ay,az},{bx,by,bz})
Matrix product
matrix({{a11,...,a1k},...,{an1,...,ank}},{{b11,...,b1p},...,{bk1,...,bkp}})
matrix({{a11,...,a1k},...,{an1,...,ank}},{b1,...,bk})
matrix({a1,...,ak},{{b11,...,b1p},...,{bk1,...,bkp}})
Tensor product
tensor({{a11,...,a1p},...,{an1,...,anp}},{{b11,...,b1q},...,{bm1,...,bmq}})
Complex and quaternion product
complex({a,b},{c,d})
quaternion({a,b,c,d},{e,f,g,h})
Geometric product and differential operator
geometric({a1,...,am},{b1,...,bm},algebra_opt)
del({f1(x1,...,xn),...,fm(x1,...,xn)},{x1,...,xn},algebra_opt)
m=2^n
algebra : cl(p,q), p+q=n
Transposition, trace, determinant
tran({{a11,...,a1p},...,{an1,...,anp}})
trace({{a11,...,a1n},...,{an1,...,ann}})
det({{a11,...,a1n},...,{an1,...,ann}})
Polynomial coefficients and solving
coef(p(x),x), solve(p(x),x,subscript_opt)
Variable substitution
subst(f(x),x,value), subst(f(x1,...,xn),{x1,...,xn},{a1,...,an})
Limit
lim(f(x),x,value,direction_opt) value : -infin, a, infin direction : <0, 0, >0
Sum, product
sum(a[i],i,n1,n2), prod(a[i],i,n1,n2)
Comparison
eq(a,b) le(a,b) ge(a,b) ne(a,b) lt(a,b) gt(a,b) ap(a,b)
Number theory
modpow(a,exponent,modulo) modinv(a,modulo) eulerphi(a) primitiveroots(a)
Groebner basis computation
groebner({p1(x1,...,xm),...,pn(x1,...,xm)},{x1,...,xm},ordering_opt,modulo_opt)
ordering : lex=lexicographic, tdl=total degree lexicographic, drl=degree reverse lexicographic, elim(k)=kth-elimination
modulo : prime integer
Expand
The expand tries to evaluate your expression. The expression below checks the output in form of strings (this they must return true):
j.expand("27^(1/3)/3") == "1"
j.expand("sqrt((2^(2^6)-1)^2)") == "18446744073709551615"
j.expand("81^(1/4)") == "3"
Converting to elementary functions
j.elementary("abs(x)") == "sqrt(x^2)"
j.elementary("sgn(x)") == "x/sqrt(x^2)"
Numeric
Reduces to numeric values:
j.numeric("pi") == "3.141592653589793"
j.numeric("1.+exp(2.)") ="8.38905609893065"
j.numeric("cos(pi)") == "-1.0"
j.numeric("Infinity") == "Infinity"
j.numeric("-1/0") == "Infinity"
j.numeric("0/0") == "NaN"
j.numeric("-1/0") == "-Infinity"
Simplify
It is often useful to rewrite an expression in term of elementary functions (log, exp, frac, sqrt, implicit roots), before simplifying it.
Below we show a few simplification examples. Such tests will return true since we compare the result of simplification with
the expectation string:
j.simplify("exp(sqrt(-1)*pi)") == "-1"
j.simplify(j.elementary("cos(x)^2+sin(x)^2")) == "1"
j.simplify("exp(sqrt(-1)*pi)") == "-1"
j.simplify(j.elementary("cos(acos(x))")) == "x"
All trigonometric functions can be reduced to "elementary" functions:
j.simplify(j.elementary("atan(tan(x))")) == "(sqrt(-1)*log(1/exp(sqrt(-1)*x)^2))/2"
j.simplify(j.elementary("cos(x)")) == "(1+exp(sqrt(-1)*x)^2)/(2*exp(sqrt(-1)*x))"
j.simplify(j.elementary("sin(x)")) == "(sqrt(-1)-sqrt(-1)*exp(sqrt(-1)*x)^2)/(2*exp(sqrt(-1)*x))"
j.simplify(j.elementary("acos(x)")) == "sqrt(-1)*log(x+sqrt(-1+x^2))"
j.simplify(j.elementary("asin(x)")) == "sqrt(-1)*log(-sqrt(-1)*x+sqrt(1-x^2))"
j.simplify(j.elementary("atan(x)")) == "(sqrt(-1)*log((sqrt(-1)+x)/(sqrt(-1)-x)))/2"
j.simplify(j.elementary("acosh(x)")) == "log(x+sqrt(-1+x^2))"
j.simplify(j.elementary("asinh(x)")) == "log(x+sqrt(1+x^2))"
j.simplify(j.elementary("atanh(x)")) == "log((1+x)/(1-x))/2"
j.simplify(j.elementary("tanh(x)")) == "-(1-exp(x)^2)/(1+exp(x)^2)"
j.simplify("x^(1/-2)") == "x^(1/-2)"
Here is yacas/wester coding. You can use "pi" instead 3.14:
j.simplify("sqrt(2*sqrt(3)+4)-(1+sqrt(3))") == "-1-sqrt(3)+sqrt(2)*sqrt(2+sqrt(3))"
j.simplify(j.elementary("cos(5*pi/6)+sin(pi/4)"))
You may represent it as:
j.simplify(j.elementary("acos(cos(x))"))
as this long expressions.
"-sqrt(-1)*log(2)+sqrt(-1)*log(sqrt((1-2*exp(sqrt(-1)*x)^2+exp(sqrt(-1)*x)^4)/exp(sqrt(-1)*x)^2) +3*exp(sqrt(-1)*x)+sqrt((1-2*exp(sqrt(-1)*x)^2 +exp(sqrt(-1)*x)^4)/exp(sqrt(-1)*x)^2)^2*exp(sqrt(-1)*x)-exp(sqrt(-1)*x)^3)"
Below is the most trivial symbolic substitutions performed by the code:
j.simplify("log(-2)") == "sqrt(-1)*pi+log(2)"
j.simplify("sqrt(0)") == "0"
j.simplify("sqrt(1)") == "1"
j.simplify("sqrt(-2)") == "sqrt(-1)*sqrt(2)"
j.simplify("sqrt(-4)") == "2*sqrt(-1)"
j.simplify("exp(0)") == "1"
j.simplify("exp(-1)") == "1/exp(1)"
j.simplify("exp(2)") == "exp(1)^2"
j.simplify("exp(-4)") == "1/exp(1)^4"
j.simplify("1/-x") == "-1/x"
Let's take a look at this example:
j.simplify("sqrt(17/12)-(1/6*sqrt(51))") == "0"
j.simplify("log(6)") == "log(2)+log(3)"
j.simplify("(x + 1)^(1/x)") == "(1+x)^(1/x)"
j.simplify("(x^2-1)/(x-1)") == "1+x"
j.simplify("abs(abs(x))" == "abs(x)"
j.simplify("abs(sgn(x))") == "1"
Examples with "conjugate"
j.simplify("conjugate(conjugate(x))" == "x"
j.simplify("conjugate(x^2)") == "conjugate(x)^2"
Substitution
j.expand("subst(1/x^2,x,a)") # output: "1/a^2"
j.expand("subst(sqrt(x+1/2)+1/sqrt(x-1/2),1/2,a)") # output 1/sqrt(-a+x)+sqrt(a+x)
j.expand("subst((a)*(b),a,c)") # "b*c"
j.expand("subst({a,b}*c,c,{x,y})") # {a, b}*{x, y}
j.expand("subst(subst(x,x,z),x,y)") # output: z
j.expand("subst((a),(a),x)") # x
j.expand("subst((a),a,x)") # x
j.expand("subst(a,(a),x)") # x
j.expand("subst(a+b,{a,b},{x,y})") # "x+y"
j.expand("subst({a,b}+{a,c},{a,b},{x,y})") # "{x, c}+{x, y}"
Substitutions can be done using differential:
j.expand("subst(x,x,d(cos(t),t))") # "d(cos(t), t)"
Substitutions can be done using matricies
j.expand(j.expand("subst(subst(matrix(A,B),A,{{a,b},{c,d}}),B,{x,y})")) # "{a*x+b*y, c*x+d*y}"
j.expand("subst(subst(A*B,A,{{a,b},{c,d}}),B,{x,y})")) # {x, y}*{{a, b},\n"+ "{c, d}}"
j.expand("subst({{a,b},{c,d}},a,A)") # {{A, b},\n"+"{c, d}}"
j.expand("subst({{a,b},{c,d}}*e,a,A)") # e*{{A, b},\n"+"{c, d}}
j.expand("subst(matrix({{a,b},{c,d}},{a,b}),{{a,b},{c,d}},A)") # "matrix(A, {a, b})"
Differentiate
j.expand("d(cos(f(x)),x)") # output "-f'(x)*sin(f(x))"
Integrate
Symbolic integration is done using "integral" operator. Use "expand" function to get the results.
Here are indefinite integrals:
j.expand("integral(1/(1+x^2),x)") == "root[0](1, 0, 4)*log(1-2*x*root[0](1, 0, 4))+root[1](1, 0, 4)*log(1-2*x*root[1](1, 0, 4))"
j.expand("integral(x+7/5,x)") == "7/5*x+1/2*x^2"
You can use any arbitrary parameter, as "a" in this example:
j.expand("integral(1/(x+a),x)") == "log(a+x)"
You can cross check the result by differentiating the output:
j.simplify(expand("d("log(a+x)",x)")) == "1/(a+x)"
Or yo can do "definite" by specifying the range:
j.expand("integral(1/x,x,1,10)")) == "log(10)"
There are more complicate examples:
j.expand("integral(tan(a+b*x),x)") == "-1/b*log(4*cos(a+b*x))"
j.simplify(j.expand("integral((x^4 + 2*x^2 + 2*x + 1) / x^3, x)")) == "-(1+4*x-x^4-4*x^2*log(x))/(2*x^2)"
Outputs =
j.toMathML("x^2*y^2"), "<math><mrow><msup><mi>x</mi><mn>2</mn></msup><msup><mi>y</mi><mn>2</mn></msup></mrow></math>"
Polynomial Solving
The engine allows polynomial solving:
j.expand("solve(c+b*x+a*x^2,x)") # answer: root[0](c, b, a)
Algebraic equation systems
One solve algebraic equation systems of any degree, with several indeterminates, by computing the Groebner bases of polynomial ideals. For example, let this system for the indeterminates x, y:
x^2 + y^2 = 4 x*y = 1
We use:
j.expand("groebner({x^2 + y^2 - 4, x*y - 1}, {x, y})")
The returned output is:
{1-4*x^2+x^4, 4*x-x^3-y}
which allows to find x, then y from x. This operation doesn't calculate the roots, it just writes the equation. For example, it wouldn't give "a = 4/5" but "5*a-4" ("= 0" implied). Groebner basis computation is explained in more details below:
Math background
I - Groebner basis computation
Here are some theoretical data gathered from [1]: An ideal V is the set of all polynomials g*f, g = any polynomial, f belonging to a family F of generator polynomials. We consider the polynomials in the indeterminates x1,x2,...,xn. Let < be an ordering over the monomials that satisfies to the following conditions: a) if a < b, then for every monomial c, we have ac < bc b) for every monomials a and b with b!=1, we have a < ab
The lexicographical ordering satisfies to these conditions. It is defined as:
x1^a1*x2^a2*...*xn^an < x1^b1*x2^b2*...*xn^bn if: aj = bj for j<i ai<bi
The maximum monomial of a polynomial with respect to the ordering < is called its head monomial.
Let G be a set of polynomials. We say that a polynomial f is reduced against G if no head monomial of a polynomial of G divides the head monomial of f. If f isn't reduced, we can subtract multiples of the elements of G from it. This process is called a reduction of f against G. It is finite, i.e. it ends up with a reduced polynomial.
A family of generators (a "basis") G of an ideal I is called a Groebner basis with respect to the ordering < if the reduction against G of every polynomial f of I always gives 0.
Theorem: A system of polynomial equations has a finite number of solutions (it is "zero dimensional") if and only if every indeterminate appears alone in the head monomial of one polynomial of the corresponding Groebner basis with respect to any ordering satisfying to the conditions a) and b) above.
When the basis is computed with respect to an "elimination" ordering (for instance, lexicographical), this gives a method to compute every solution of the system by finding one indeterminate after the other. (This is one motivation for computing a Groebner basis).
To compute a Groebner basis, we can use Buchberger's algorithm: We define the s-polynomial of two polynomials f and g as:
S(f,g)=(scm(hm(f),hm(g))/hm(f))*f - (scm(hm(f),hm(g))/hm(g))*g
where hm(f) (resp hm(g)) stands for the head monomial of f (resp g).
Theorem: A basis G is a Groebner basis if and only if for every pair (f,g) of polynomials of G, the reduction of S(f,g) against G gives 0.
If some S(f,g) doesn't reduce to 0, we can add it to G (since it's a linear combination of polynomials of G, the ideal generated by G remains the same). S(f,g) then reduces to 0 against the new G. We can repeat this process for every S(f,g). Buchberger has proved that this process is finite, i.e. it ends up with a Groebner basis.
At each step, we must make the choice of a pair (f,g). This is where some optimization can take place. We can eliminate some possibilities using Buchberger's third criterion:
If there is an element h of G such that the head monomial of h divides the scm of the head monomials of two polynomials f and g, and if we already have considered S(f,h) and S(h,g), then we don't have to consider S(f,g) because it reduces to 0.
Even after applying this criterion, several possibilities may remain. According to Buchberger, the choice of a pair such that the scm of the head monomials is minimal in the current ordering is good. As an improvement, we chose the pair with the minimal "sugar" (see [3]).
II - Symbolic integration
The rational functions integration algorithm is taken from [4].
III - Number theory
There are good explanations in wikipedia which I used when implementing Euler's phi function and primitive roots.
References
[1] Davenport J., Siret Y., Tournier E., Calcul formel - systemes et algorithmes de manipulations algebriques, MASSON, 1993 ( http:staff.bath.ac.uk/masjhd/masternew.pdf )
[2] JC Faugere ( http:calfor.lip6.fr/~jcf/ )
[3] Alessandro Giovini, Teo Mora, Gianfranco Niesi, Lorenzo Robbiano and Carlo Traverso, "One sugar cube, please" or selection strategies in the Buchberger algorithm, Proceedings of the 1991 international symposium on Symbolic and algebraic computation July 15 - 17, 1991, Bonn West Germany
[4] Manuel Bronstein, Symbolic integration tutorial, ISSAC'98, Rostock, August 12, 1998 ( http:www-sop.inria.fr/cafe/Manuel.Bronstein/publications/issac98.ps.gz )
[5] Thom Mulders, A note on subresultants and a correction to the lazard/rioboo/trager formula in rational function integration, Journal of Symbolic Computation, 24(1):45-50, 1997
[6] F Mancini ( http:www.geocities.com/famancin/groebner_applet.html )
[7] Singular ( http:www.singular.uni-kl.de )
[8] David Hestenes, Oersted Medal Lecture 2002: Reforming the Mathematical Language of Physics ( http:modelingnts.la.asu.edu/pdf/OerstedMedalLecture.pdf )
[9] The Geobucket Data Structure for Polynomials, T Yan - JSC, 1998 ( http:historical.ncstrl.org/tr/ps/cornellcs/TR96-1607.ps )
