容易计算
$$s_n=900\cdot\frac{1-r^n}{1-r}-3(\frac{1-r^{n+1}}{1-r})'=900\cdot\frac{1-r^n}{1-r}-3\cdot\frac{-(n+1)r^n(1-r)+(1-r^{n+1})}{(1-r)^2},$$
化简得
$$900(1-r^n)(1-r)+3(n+1)r^n(1-r)-3(1-r^{n+1})-s_n(1-r)^2=0.$$
即
$$900(1-r^n)+3(n+1)r^n-3(1+r+\ldots+r^n)-s_n(1-r)=0.$$
代入$n=5000$和$s_n=-6\times 10^{11}$得
$$14100r^{5000}-3(r^{4999}+\ldots+r^2)-600000000003r+600000000897=0.$$
将上面的等式左侧记作$f(r)$。
用geogebra的图像计算器可以画出f的图像，并将根计算到小数点后10位，即$1.0023221086$。
接下来以此为初值，用牛顿法
$$x_{k+1}=x_k-\frac{f(x)}{f'(x)},$$
迭代计算，当相邻两次迭代之间的差异小于$10^{-13}$时停止。
由于$f$的次数较高，此处用秦九韶法计算多项式的值，即
$$a_nx^n+\ldots+a_0=(\ldots((a_nx+a_{n-1})x+a_{n-2})x+\ldots+a_1)x+a_0.$$
等式右侧的计算量是线性的，可以用一个栈处理。
注：本题中涉及的秦九韶法、牛顿法一般在《数值分析》等课程中讲授。第一个公式（即$s_n$的公式）中的导数部分是我国《高等数学》课程的标准内容。
import decimal
decimal.getcontext().prec=15
def p(a,x):
b=a[:]
while len(b)>1:
u,v=b.pop(),b.pop()
b.append(u*x+v)
return b[0]
f=[600000000897,-600000000003]+[-3]*4998+[14100]
df=[-600000000003]+[-3*i for i in range(2,5000)]+[14100*5000]
epsilon=decimal.Decimal(0.0000000000001)
x0=decimal.Decimal(1.0023221086)
x1=x0-p(f,x0)/p(df,x0)
while abs(x1-x0)>=epsilon:
x0=x1
x1=x0-p(f,x0)/p(df,x0)
decimal.getcontext().prec=13
print decimal.Decimal(1.0)*x1
|
