Ask Your Question

Solving a power series equation

asked 2018-06-29 12:52:53 -0500

Daniel L gravatar image

I'd like to solve the equation f(t) == 1 where f is a power series:

a1 = float(4 + sqrt(8)) 
b1 = float(3 + sqrt(8)) 
c1 = b1
an = 4*n*(n-1)+a1*n
an1 = 4*n*(n+1) + a1*(n+1)
anr = b1 + 2*n*(4*n+a1)
bn = 1+2*n*(n-1)+n*(b1-1)
bn1 = 1+2*n*(n+1)+(n+1)*(b1-1)
bnr = c1 + 2*n*(2*n+b1-1)

def f(t):
return (1+b1)^-t + 3/2*a1^-t + 1/2*(2*b1)^-t + b1^-t  + 2*sum(an^-t +an1^-t + anr^-t + 2*bn^-t +2*bn1^-t + bnr^-t ,n,1,1000).n()

I've tried simply plugging in values for t and I estimate that if f(t) = 1 then t is approximately 1.51. I'd like to possibly do something like

solve(f(t) == 1,t)

but this appears to lock sage up and I have to interrupt the process. The true f function is actually an infinite series, but I am truncating it to the first 1000 terms. I know there are methods of solving series equations such as the method of regula falsi. Does Sage have anything like this? Thanks!

edit retag flag offensive close merge delete


If an answer solves your question, please accept it by clicking the "accept" button (the one with a check mark, at the top left of the answer).

This will mark the question as solved in the Ask Sage's list of questions, as well as in lists of questions related to a particular query or keyword.

slelievre gravatar imageslelievre ( 2018-07-02 03:11:25 -0500 )edit

1 answer

Sort by ยป oldest newest most voted

answered 2018-06-29 13:29:51 -0500

updated 2018-07-01 05:42:17 -0500

What your code for f(t) does is first take the sum symbolically up to 1000 terms, and then evaluate numerically.

What you want to do instead is calculate the terms numerically and sum them.

For example, like this:

def f(t):
    return (1+b1)^-t + 3/2*a1^-t + 1/2*(2*b1)^-t + b1^-t  + 2*sum((an^-t +an1^-t + anr^-t + 2*bn^-t +2*bn1^-t + bnr^-t).subs(n=k).n() for k in xrange(1,1000))

The method solve is used to find exact solutions. For approximate solutions you can use e.g. find_root:

find_root(lambda t: f(t) - 1, 1.48, 1.52)

It yields 1.5111899890344738, and f(1.5111899890344738) is approximately 1.00000000000000.

edit flag offensive delete link more


When I run this code I get the following traceback error: ValueError: rtol too small (4.5e-16 < 8.88178e-16) I'm not sure what this means or why it's happening.

Daniel L gravatar imageDaniel L ( 2018-06-30 15:58:02 -0500 )edit

That seems to be (a bug) in some version(s) of Sage. You can work around it by specifying the rtol keyword argument to find_root explicitly, or using the partial method in the thread I just linked, or (I would hope) by upgrading to the latest version of SageMath.

rburing gravatar imagerburing ( 2018-06-30 16:25:18 -0500 )edit

The workaround on the above thread works well. Thanks for the help!

Daniel L gravatar imageDaniel L ( 2018-06-30 17:13:48 -0500 )edit

You're welcome. You can mark this answer as correct.

rburing gravatar imagerburing ( 2018-07-01 05:46:26 -0500 )edit

Your Answer

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

Add Answer

Question Tools

1 follower


Asked: 2018-06-29 12:52:53 -0500

Seen: 98 times

Last updated: Jul 01