DMelt:Symbolic/3 JScl

From HandWiki
Revision as of 09:55, 26 December 2016 by imported>Jworkorg
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Member

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 )