How to iterate over symbolic coefficients

asked 2019-10-18 10:29:04 +0100

Edgar Brown
79 ●3 ●3 ●12

updated 2019-10-28 01:01:56 +0100

I am a relative newbie to Sage but not to algebraic manipulation in general.

I have a family of expressions with non-linear coefficients which I am not sure how best to represent. In the end it has to be an (infinite) ring of polynomials representing a Volterra space, but I have the same terms inside non-linear functions (which I need to evaluate to a known value) and as the polynomial variables.

Here is an example, where simple symbolic manipulation doesn't work:

f_m = a*x^b
ex = expand(f_m.subs({b : b_m + (b_s-b_m)^3})).derivative(b_m).full_simplify()
ex.subs({b_m + (b_s-b_m)^3 : b_m }).show()

Producing:

  -(3*a*b_m^2 - 6*a*b_m*b_s + 3*a*b_s^2 - a)*x^(-b_m^3 + 3*b_m^2*b_s - 3*b_m*b_s^2 + b_s^3 + b_m)*log(x)

When what I want is:

  -(3*a*b_m^2 - 6*a*b_m*b_s + 3*a*b_s^2 - a)*x^(-b_m)*log(x)

How can I iterate over non-linear elements (exponent and the arguments to any non-monomial terms) carrying out direct substitutions, while leaving monomials alone. I want the result to remain as a polynomial ring on b_s

Clarification: The intent is to use the code for a generic non-linear construct. It's intended to handle any general non-linear expression in L2, but these are manipulated to belong to a class of Volterra polynomials with non-linear coefficients. like the log(x) and x^bm on the above example. Those are the coefficients of the monomials.

Idealy whatever process I follow should be able to identify the non-monomial coefficients in the expression and remove the monomial terms from them by doing the substitution b_s = b_m. Basically:

    x.subs(b_s = b_m) for x not a monomial term in expression.

Which would make the resulting expression strictly a polynomial ring on b_s with non-linear coefficients on x and b_m.

edit retag flag offensive close merge delete

Comments

Are you using the "symbolic ring" or the more specific polynomial rings? It's unclear from your code, though the use of full_simplify indicates you are just using var() to create new variables - making a polynomial ring in b_s might help, I'm not sure.

kcrisman ( 2019-10-18 13:47:27 +0100 )edit

I don't have experience with either, but it's obvious to me that a polynomial ring is what I want and a symbolic ring might provide me with the proper iteration tools. When I tried the polynomial ring of polynomial coefficients I kept getting errors due to the same variable being used on the exponents and the monomials. I have no idea about the symbolic ring or how to use it, but given that this is just a minor step in a long process doing things by hand is not an option.

Edgar Brown ( 2019-10-18 15:09:41 +0100 )edit

Are b_m and b_s supposed to remain unknown, or will they be given actual integer values ?

tmonteil ( 2019-10-18 15:28:51 +0100 )edit

@tmontell They are supposed to remain unknown until I can coerce the final output of a long set of operations into a numeric framework.

Edgar Brown ( 2019-10-18 15:41:26 +0100 )edit

add a comment

3 Answers

Sort by » oldest newest most voted

answered 2019-10-19 00:17:25 +0100

tmonteil
27423 ●31 ●204 ●517 http://wiki.sagemath.o...

Okay, i am not sure to understand your question in whole generality. However, here is a solution that could be generalized and adapted to your needs. Let me assume that you want to rewrite the exponents of your expression and that those exponents are polynomials. Working on polynomials has the good property that instead of doing some pattern matching as in the symbolic ring, you can use quotient rings.

For example, here is what you can do in your previous example:

sage: R.<b_m,b_s> = PolynomialRing(QQ)
sage: R
Multivariate Polynomial Ring in b_m, b_s over Rational Field
sage: S = R.quotient(b_m + (b_s-b_m)^3 - b_m)
sage: S
Quotient of Multivariate Polynomial Ring in b_m, b_s over Rational Field by the ideal (-b_m^3 + 3*b_m^2*b_s - 3*b_m*b_s^2 + b_s^3)
sage: p = -b_m^3 + 3*b_m^2*b_s - 3*b_m*b_s^2 + b_s^3 + b_m
sage: p
-b_m^3 + 3*b_m^2*b_s - 3*b_m*b_s^2 + b_s^3 + b_m
sage: S(p)
b_mbar
sage: lift(S(p))
b_m

Now, what you have to do is to walk along the expression and apply this recipe to all exponents of the powers that are encountered. A power is a kind of arithmetic operation, so you will have to override the arithmetic method of ExpressionTreeWalker. Here is how to:

from sage.symbolic.expression_conversions import ExpressionTreeWalker
class RewritePolynomialPowers(ExpressionTreeWalker):
    def __init__(self, ring):
        self.ring = ring
    def arithmetic(self, ex, op):
        if op is operator.pow:
            base, exponent = ex.operands()
            exponent = SR(str(lift(self.ring(exponent))))
            return base ** exponent
        else:
            return reduce(op, map(self, ex.operands()))

Let us try with your example:

sage: var('a,x,b,b_s,b_m')
(a, x, b, b_s, b_m)
sage: f_m = a*x^b
sage: expr = expand(f_m.subs({b : b_m + (b_s-b_m)^3})).derivative(b_m).full_simplify()
sage: RewritePolynomialPowers(S)(expr)
-(3*a*b_m^2 - 6*a*b_m*b_s + 3*a*b_s^2 - a)*x^b_m*log(x)

Please tell if this is sufficient for your general case. Anyway, you see the TreeWalker strategy, you can also replace other kind of operators (composition, derivation, relation,...), see the (not documented) https://doc.sagemath.org/html/en/refe...

edit flag offensive delete link

Comments

Yikes, thanks!! It will take me a while to digest that, I was hoping for some simpler pattern matching.

Could this be easily extended to other non-linear functions such as log(x), sin(x), etc.? Or particular cases would have to be treated one by one (e.g., by iterating through a list of all nonlinear functions). I don't know enough to see if there are edge cases, but I will definitively look into TreeWalker.

Edgar Brown ( 2019-10-19 04:02:04 +0100 )edit

VERY nice.

This TreeWalker should have some mention in the doc, and a docstring... It has potentially a zillion uses in "RealWorld"(TM) problems.

Emmanuel Charpentier ( 2019-10-19 10:45:12 +0100 )edit

Indeed ! I was thinking of a tutorial instead of just doc, since the way it works is by overriding a default class, not by calling a function. This is a standard pattern in some Python programs, but not in Sage. In particular, i had to copy a part of the source code to insert my hook, so i guess it is not very handy. Perhaps should we propose a better user interface.

tmonteil ( 2019-10-19 11:36:14 +0100 )edit

I might have found a much simpler and general way to solve my specific issue that just involves a couple of function manipulations, I'll try to implement it and post it as soon as I get a chance by the end of this coming week.

Edgar Brown ( 2019-10-20 05:26:09 +0100 )edit

I believe your answer provides the path to the only reliable choice. I think I have to walk the tree of what I assume to be mostly a polynomial looking at sums and products of elements and expanding the substitution in all of those cases, but stopping at anything else such as division and exponentiation and removing the substitution on those subtrees. Although conceptually it works, I can't see the path for doing that in code. Your example enumerates the tree backwards from this, which means that I would have to account for all (unknown) cases, instead of known ones such as sums and products.

I still get the feeling that a properly defined symbolic ring would be able to do this, but I don't see how.

Edgar Brown ( 2019-10-28 01:01:09 +0100 )edit

see more comments

answered 2019-10-25 14:33:56 +0100

Edgar Brown
79 ●3 ●3 ●12

updated 2019-10-28 00:53:10 +0100

Leaving this answer as a warning. I decided to go in the direction of @tmontell 's answer, although I can see multiple ways to get the functions to do what I want in specific instances, there are too many corner cases that will cause problems in the future.

Although @tmonteil answer walks the structure, I believe this is simpler to implement and understand. This answer does what I need, at least for the cases that I am handling by hand (instead of programmatically):

WARNING: There is an issue where if I simplify before applying the final substitution these do not work, I can't see a reason for that. This question might illuminate the reason.

Definitions:

  var('a b x b_m b_s',domain='real')
  faux = function('faux', latex_name='f_{aux}')(b_m) # Auxiliary function
  R = PolynomialRing(SR,names='a,b_s')  # Polynomial ring over the symbolic ring
  f_m = a*x^b     # target function

Note that the definitions of R is critical, a & b_s are not the same as the variables being used (see answer: map functions into polynomial coefficients).

The actual operations:

 exd = expand(f_m.subs(b = faux)).derivative(b_m)  # substitute the auxiliary function
 exd = exd.derivative(b_m)      # calculate the derivative
 exd = exd.subs({diff(faux,b_m) : 1 - 3*(b_s-b_m)^2})  # replace by known derivative 
 exd = exd.subs({faux : b_m}).expand()    # replace by function evaluation

This set of operations do what I need, now all that is left is to transform exd into a proper polynomial so I can manipulate its coefficients.

 exp = exd.polynomial(ring=R)         # It's now represented as a polynomial
 dict(zip(exp.monomials(), exp.coefficients()))  # Dictionary with the desired terms

Which results in:

 {a: -3*b_m^2*x^b_m*log(x) + x^b_m*log(x),
  a*b_s: 6*b_m*x^b_m*log(x),
  a*b_s^2: -3*x^b_m*log(x)}

Some cleanup might be needed, and there are some issues (e.g., PolynomialRing does not support "show()" for some reason), and in some cases the substitution breaks for no apparent reason (a bug?). But I think it does what is needed. Any suggestions on how to improve this?

edit flag offensive delete link

add a comment

answered 2019-10-18 13:48:53 +0100

Emmanuel Charpentier
7764 ●6 ●50 ●153

updated 2019-10-18 13:51:03 +0100

I probably do not understand your problem: your desired expression can be quite simply obtained by abstaining to expand and full_simplify the original expression:

sage: f_m.subs({b : b_m + (b_s-b_m)^3}).diff(b_m)
-(3*(b_m - b_s)^2 - 1)*a*x^(-(b_m - b_s)^3 + b_m)*log(x)

On the other hand, you may expand the left hand of your second substitution :

sage: ex
-(3*a*b_m^2 - 6*a*b_m*b_s + 3*a*b_s^2 - a)*x^(-b_m^3 + 3*b_m^2*b_s - 3*b_m*b_s^2 + b_s^3 + b_m)*log(x)
sage: ex.subs({(b_m + (b_s-b_m)^3).expand(): b_m})
-(3*a*b_m^2 - 6*a*b_m*b_s + 3*a*b_s^2 - a)*x^b_m*log(x)

edit flag offensive delete link

Comments

I know that with playing around with the expressions I can get Sage to do it, but this is just a minor step in a complicated process, so I need a reliable and automatable method that goes beyond simple pattern matching. Iterating though the non-monomial elements seems to be the ticket, but I don't know how to do that.

Edgar Brown ( 2019-10-18 15:08:10 +0100 )edit

Could you please provide a more general form of expressions you will have to handle ?

tmonteil ( 2019-10-18 15:41:05 +0100 )edit

@tmontel I added some clarification to the question. I hope that covers it.

Edgar Brown ( 2019-10-18 22:22:27 +0100 )edit

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

How to iterate over symbolic coefficients

Comments

3 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

How to iterate over symbolic coefficients edit

Comments

3 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

How to iterate over symbolic coefficients