Taylor expansion in SageMath returns error

Taylor_expansion

asked 2024-08-13 18:27:57 +0200

AlexCg
1 ●1

updated 2024-08-14 10:36:14 +0200

Hi, I need to find the Taylor coefficient of order 7 in the variable "x4" of the following function:

var('x1 x2 x3 x4 w1 w2 w3 w4 w5 w5 w6');

g= 1/4*((w1^2*w3^2*x1^4*x2^6*x3^8 + 2*(w1^2*w3^2*x1^4*x2^6 + w1^2*w3*x1^3*x2^5)*x3^7 + (w1^2*w3^2*x1^4*x2^6 + 2*w1^2*w3*x1^3*x2^5 + 2*w1^2*x1^2*x2^4)*x3^6)*x4^6 + 2*((w1^2*w3^2*x1^4*x2^6 + (w1^2*w3^2*x1^4 + (2*w1^2*w3 + w1*w3^2)*x1^3)*x2^5)*x3^7 + 2*(w1^2*w3^2*x1^4*x2^6 + (w1^2*w3^2*x1^4 + (3*w1^2*w3 + w1*w3^2)*x1^3)*x2^5 + (w1^2*w3*x1^3 + (w1^2 + w1*w3)*x1^2)*x2^4)*x3^6 + (w1^2*w3^2*x1^4*x2^6 + (w1^2*w3^2*x1^4 + (4*w1^2*w3 + w1*w3^2)*x1^3)*x2^5 + 2*(w1^2*w3*x1^3 + (2*w1^2 + w1*w3)*x1^2)*x2^4 + 2*(w1^2*x1^2 + w1*x1)*x2^3)*x3^5)*x4^5 + ((w1^2*w3^2*x1^4*x2^6 + 2*(w1^2*w3^2*x1^4 + (3*w1^2*w3 + 2*w1*w3^2)*x1^3)*x2^5 + (w1^2*w3^2*x1^4 + 2*(3*w1^2*w3 + 2*w1*w3^2)*x1^3 + 2*(3*w1^2 + 6*w1*w3 + w3^2)*x1^2)*x2^4)*x3^6 + 2*(w1^2*w3^2*x1^4*x2^6 + (2*w1^2*w3^2*x1^4 + (7*w1^2*w3 + 4*w1*w3^2)*x1^3)*x2^5 + (w1^2*w3^2*x1^4 + 4*(2*w1^2*w3 + w1*w3^2)*x1^3 + (9*w1^2 + 16*w1*w3 + 2*w3^2)*x1^2)*x2^4 + (w1^2*w3*x1^3 + (3*w1^2 + 4*w1*w3)*x1^2 + 2*(3*w1 + w3)*x1)*x2^3)*x3^5 + (w1^2*w3^2*x1^4*x2^6 + 2*(w1^2*w3^2*x1^4 + 2*(2*w1^2*w3 + w1*w3^2)*x1^3)*x2^5 + (w1^2*w3^2*x1^4 + 2*(5*w1^2*w3 + 2*w1*w3^2)*x1^3 + 2*(7*w1^2 + 10*w1*w3 + w3^2)*x1^2)*x2^4 + 2*(w1^2*w3*x1^3 + (5*w1^2 + 4*w1*w3)*x1^2 + 2*(5*w1 + w3)*x1)*x2^3 + 2*(w1^2*x1^2 + 4*w1*x1 + 2)*x2^2)*x3^4)*x4^4 + 2*(((w1^2*w3 + w1*w3^2)*x1^3*x2^5 + (2*(w1^2*w3 + w1*w3^2)*x1^3 + (3*w1^2 + 10*w1*w3 + 2*w3^2)*x1^2)*x2^4 + ((w1^2*w3 + w1*w3^2)*x1^3 + (3*w1^2 + 10*w1*w3 + 2*w3^2)*x1^2 + 4*(3*w1 + 2*w3)*x1)*x2^3)*x3^5 + (2*(w1^2*w3 + w1*w3^2)*x1^3*x2^5 + (4*(w1^2*w3 + w1*w3^2)*x1^3 + (7*w1^2 + 22*w1*w3 + 4*w3^2)*x1^2)*x2^4 + 2*((w1^2*w3 + w1*w3^2)*x1^3 + 2*(2*w1^2 + 6*w1*w3 + w3^2)*x1^2 + (17*w1 + 10*w3)*x1)*x2^3 + ((w1^2 + 2*w1*w3)*x1^2 + 2*(5*w1 + 2*w3)*x1 + 8)*x2^2)*x3^4 + ((w1^2*w3 + w1*w3^2)*x1^3*x2^5 + 2*((w1^2*w3 + w1*w3^2)*x1^3 + (2*w1^2 + 6*w1*w3 + w3^2)*x1^2)*x2^4 + ((w1^2*w3 + w1*w3^2)*x1^3 + (5*w1^2 + 14*w1*w3 + 2*w3^2)*x1^2 + 12*(2*w1 + w3)*x1)*x2^3 + ((w1^2 + 2*w1*w3)*x1^2 + 2*(7*w1 + 2*w3)*x1 + 12)*x2^2 + 2*(w1*x1 + 2)*x2)*x3^3)*x4^3 + 24*(x2^2 + 2*x2 + 1)*x3^2 + 2*(((w1^2 + 4*w1*w3 + w3^2)*x1^2*x2^4 + 2*((w1^2 + 4*w1*w3 + w3^2)*x1^2 + (8*w1 + 7*w3)*x1)*x2^3 + ((w1^2 + 4*w1*w3 + w3^2)*x1^2 + 2*(8*w1 + 7*w3)*x1 + 20)*x2^2)*x3^4 + 2*((w1^2 + 4*w1*w3 + w3^2)*x1^2*x2^4 + (2*(w1^2 + 4*w1*w3 + w3^2)*x1^2 + 3*(6*w1 + 5*w3)*x1)*x2^3 + ((w1^2 + 4*w1*w3 + w3^2)*x1^2 + 4*(5*w1 + 4*w3)*x1 + 27)*x2^2 + ((2*w1 + w3)*x1 + 7)*x2)*x3^3 + ((w1^2 + 4*w1*w3 + w3^2)*x1^2*x2^4 + 2*((w1^2 + 4*w1*w3 + w3^2)*x1^2 + 2*(5*w1 + 4*w3)*x1)*x2^3 + ((w1^2 + 4*w1*w3 + w3^2)*x1^2 + 6*(4*w1 + 3*w3)*x1 + 36)*x2^2 + 2*((2*w1 + w3)*x1 + 9)*x2 + 2)*x3^2)*x4^2 + 24*x2^2 + 48*(x2^2 + 2*x2 + 1)*x3 + 12*(((w1 + w3)*x1*x2^3 + (2*(w1 + w3)*x1 + 5)*x2^2 + ((w1 + w3)*x1 + 5)*x2)*x3^3 + (2*(w1 + w3)*x1*x2^3 + (4*(w1 + w3)*x1 + 11)*x2^2 + 2*((w1 + w3)*x1 + 6)*x2 + 1)*x3^2 + ((w1 + w3)*x1*x2^3 + 2*((w1 + w3)*x1 + 3)*x2^2 + ((w1 + w3)*x1 + 7)*x2 + 1)*x3)*x4 + 48*x2 + 24)*e^(-w2*x1*x2*x3*x4)/((x2^3*x3^6*e^(w5*x1*x2) + 3*x2^3*x3^5*e^(w5*x1*x2) + 3*x2^3*x3^4*e^(w5*x1*x2) + x2^3*x3^3*e^(w5*x1*x2))*x4^3*e^(w6*x1*x2*x3) + 3*((x2^3 + x2^2)*x3^5*e^(w5*x1*x2) + 3*(x2^3 + x2^2)*x3^4*e^(w5*x1*x2) + 3*(x2^3 + x2^2)*x3^3*e^(w5*x1*x2) + (x2^3 + x2^2)*x3^2*e^(w5*x1*x2))*x4^2*e^(w6*x1*x2*x3) + 3*((x2^3 + 2*x2^2 + x2)*x3^4*e^(w5*x1*x2) + 3*(x2^3 + 2*x2^2 + x2)*x3^3*e^(w5*x1*x2) + 3*(x2^3 + 2*x2^2 + x2)*x3^2*e^(w5*x1*x2) + (x2^3 + 2*x2^2 + x2)*x3*e^(w5*x1*x2))*x4*e^(w6*x1*x2*x3) + ((x2^3 + 3*x2^2 + 3*x2 + 1)*x3^3*e^(w5*x1*x2) + 3*(x2^3 + 3*x2^2 + 3*x2 + 1)*x3^2*e^(w5*x1*x2) + 3*(x2^3 + 3*x2^2 + 3*x2 + 1)*x3*e^(w5*x1*x2) + (x2^3 + 3*x2^2 + 3*x2 + 1)*e^(w5*x1*x2))*e^(w6*x1*x2*x3))

For this we used the command:

taylor(g,x4,0,7)

Then, this return:

RuntimeError: ECL says: THROW: The catch RAT-ERR is undefined.
During handling of the above exception, another exception occurred:
......
TypeError: ECL says: THROW: The catch RAT-ERR is undefined.

However, if I compute taylor(g,x4,0,3) returns its expansion without problems, but in order > 3, returns error.

Can anybody say me why don't returns its taylor expansion in order beyond 3?

edit retag flag offensive close merge delete

Comments

Edited for legibility.

If you want to insert code in your question, indent it by four spaces (or use Ctrl-K or the "101010" button).

Emmanuel Charpentier ( 2024-08-14 10:37:55 +0200 )edit

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2024-08-14 21:25:31 +0200

Max Alekseyev

6342 ●7 ●44 ●149

updated 2024-08-21 17:18:09 +0200

You can define all variables but x4 as symbolic and then define a power series ring in variable x4 over the symbolic ring with the default precision 8 (or anything at least one more than the power of the required coefficient), the coefficient of x4^7 is then simply g[7]:

var('x1 x2 x3 w1 w2 w3 w4 w5 w6')
R.<x4> = PowerSeriesRing(SR,default_prec=8)
g = ... # copy your expression here
print( g[7] )

ADDED. From discussion with @Emmanuel Charpentier in the comments: in the setup of the code given in the question, the following addition will do the job:

R1.<p4> = PowerSeriesRing(SR,default_prec=8)
print( eval(preparse(str(g).replace('x4','p4')))[7] )

It does not require re-defining x4 as a power-series variable, but rather introduces a new variable p4 for that purpose.

edit flag offensive delete link

Comments

@Max Alekseyev : I'm afraid not, at least not directly. Exploring, using the original SR definition of g :

sage: R1.<p4>=PowerSeriesRing(SR, default_prec=8)
sage: Rg=R1(str(g).replace("x4", "p4"))
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)

[ Snip... ]

AttributeError: 'sage.symbolic.expression.Expression' object has no attribute '_SageObject__custom_name'

During handling of the above exception, another exception occurred:

[ Re-snip... ]

TypeError: -w2*x1*x2*x3*p4 is not a constant polynomial

However,

sage: %time Rg=eval(preparse("R1(%s)"%str(g).replace("x4", "p4")))
CPU times: user 1.24 s, sys: 40 µs, total: 1.24 s

Emmanuel Charpentier ( 2024-08-16 18:51:12 +0200 )edit

BTW, compare :

sage: %time Rg=eval("R1(%s)"%str(g).replace("x4", "p4").replace("^","**"))
CPU times: user 40.7 s, sys: 124 ms, total: 40.8 s
Wall time: 36.7 s

with

sage: %time Rg=eval(preparse("R1(%s)"%str(g).replace("x4", "p4")))
CPU times: user 1.24 s, sys: 40 µs, total: 1.24 s
Wall time: 1.24 s

Go figure...

Emmanuel Charpentier ( 2024-08-16 18:52:31 +0200 )edit

@Emmanuel Charpentier: I'm not sure what is your concern. It works fine at Sagecell in the way I described in my answer: https://sagecell.sagemath.org/?q=radupt

Max Alekseyev ( 2024-08-16 22:30:13 +0200 )edit

@Max Alekseyev :

I'm not sure what is your concern.

Your solution uses edited code to create a new object. I tried to generate this object from the original object, and stumbled on two coercien|conversion annoying booboos.

It works fine at Sagecell

I didn't question this. I even underscored that your solution gives an expression easier to handle than my series() solution.

Which, BTW, may hilghight a problem... I may further explore this.

Emmanuel Charpentier ( 2024-08-17 10:01:38 +0200 )edit

I see. Btw, R1() is not needed. This gives the answer:

R1.<p4> = PowerSeriesRing(SR,default_prec=8)
eval(preparse(str(g).replace('x4','p4')))[7]

Max Alekseyev ( 2024-08-18 19:51:46 +0200 )edit

see more comments

answered 2024-08-14 13:54:15 +0200

Emmanuel Charpentier
7532 ●6 ●50 ●146

updated 2024-08-22 13:09:43 +0200

taylor uses Maxima's algorithm, which may be buggy. You may explore Sage's github issues to see if this has been already reported, and possibly report it, possibly after querying sage-support.

A workaround is to use the series method of symbolic expressions, which invisibly use Sage's power series rings. Your coefficient can be computed as :

# Note that you need to ask for an 8th order series to get a significant x4^7 term !
c7=g.series(x4==0, 8).truncate().coefficient(x4, 7)
CPU times: user 780 ms, sys: 0 ns, total: 780 ms
Wall time: 779 ms

which I won't print, because it's a sum of 228 (large) terms :

sage: c7.operator()
<function add_vararg at 0x7efd2aa99620>
sage: len(c7.operands())
228

Another workaround is to use other interpreter's relevant series. For example :

# Giac
sage: %time Gc7=g._giac_().series(x4, 0, 8)._sage_().coefficient(x4, 7)
CPU times: user 274 ms, sys: 0 ns, total: 274 ms
Wall time: 492 ms
# Mathematica (here, the gratis-but-not-free Wolfram Engine)
sage: %time Mc7=g._mathematica_().Series((x4, 0, 8)).Normal().sage().coefficient(x4, 7)
CPU times: user 72 ms, sys: 632 µs, total: 72.6 ms
Wall time: 98.4 ms

Sympy's series method seems to have a hard time with this expression : it does not return in reasonable time.

Formally checking the equality of these expressions seems hard : Sage doesn't return in reasonable time. Numerically checking them may be easier but do not give a proof...

EDIT : The expressions returned by Giac and Mathematica can be proved equal :

sage: %time (Gc7-Mc7).is_zero()
CPU times: user 25.6 s, sys: 129 ms, total: 25.7 s
Wall time: 23.9 s
True

Similarly, the quantities returned by Mathematica and Sage can be proved equal by using Mathematica :

sage: %time (c7-Mc7)._mathematica_().FullSimplify().sage().is_zero()
CPU times: user 16.1 s, sys: 12 ms, total: 16.1 s
Wall time: 49.1 s
True

However, %time (c7-Mc7).is_zero() does not return in a "reasonable" time (a few minutes).

EDIT 2 : The solution proposed by @Max Alekseyev gives an expression much easier to handle by Sage's default algorithms :

sage: R1.<p4>=PowerSeriesRing(SR, default_prec=8)
sage: %time Rg=eval(preparse("R1(%s)"%str(g).replace("x4", "p4"))) # See comments below Max's anwer
sage: %time Pc7=Rg[7]
CPU times: user 100 µs, sys: 0 ns, total: 100 µs
Wall time: 113 µs
sage: %time (Pc7-Mc7).full_simplify()
CPU times: user 8.21 s, sys: 16.1 ms, total: 8.23 s
Wall time: 7.23 s
0
sage: %time (Pc7-Gc7).full_simplify()
CPU times: user 18.5 s, sys: 63 ms, total: 18.5 s
Wall time: 16.9 s
0

EDIT 2 : The Maxima problem seems attributable to Sage's interface to Maxima, since the code runs without triggering an error in a standalone Maxima session. See this sage-support thread

HTH,

edit flag offensive delete link

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

Taylor expansion in SageMath returns error

Comments

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Taylor expansion in SageMath returns error edit

Comments

2 Answers

Comments

Your Answer

Question Tools

Stats

Related questions

Taylor expansion in SageMath returns error