用python算微积分及牛顿迭代求解高阶方程
前段时间，工作中遇到了一个难题，网上查了很多资料也没有明确的解决方案，后来自己想到一种解决方法，决定分享出来。
下面先说下问题：如何计算irr(内部收益率)
我们先来看下百度百科对这1名词的解释和计算方式的描述：
看完上面的解释，给了一大堆公式，还是没有告诉怎么求，查看其它资料也大部分言尽于此，固执的我看了更多的名词解释，终于明白了这个是怎么来的，大概就是下面这个公式：
上面公式中的y和a1到an都是已知量，x就是要求的irr。这个公式怎么来的，或者说在生活中的具体例子是什么呢，可以这样理解：
小明往银行存钱，假设月利率是1%（只为举例，实际生活中肯定是低于此值的），1月到12月每月存入100块，问到年底小明能拿回多少钱，肯定是像下图这样算。
那么问题再升级一下，如果小明每月存入的金额不固定了，用an表示每月存入的金额，那结果怎么算呢？那其实还是很简单，把上面公式的100改成a1到a12就好了，如下公式：
那么问题再次升级，如果告诉你小明年底能拿到1000元，每期存的金额分别是a1到a12，假设利率是稳定不变的，此时问你利率是多少，你该怎么算呢？首先肯定还是列出下面的公式：
公式是有了，但是这是一个高阶方程，该如何求解呢，纵观数学史上，方程一旦高于5次是没有求根公式的，但是伟大的先贤们还是找到了一些解决方案，就是通过多次迭代之后，近似的找到方程的解。今天要介绍的就是几种迭代方案中的一种——牛顿迭代，用于求解f(x)=0的解：
下面我们来看看牛顿迭代的原理：
首先我们根据泰勒公式，我们知道f(x)可由其邻域x0的n阶导数的表达式近似表示，公式如下：
因为公式中后面项的分母越来越大，其对f(x)值的影响就越来越小，所以我们取线性部分，再通过迭代的方式可达到近似求解的目的：
1、取泰勒展开的线性部分
2、解线性部分的方程
3、对第2部不断进行迭代，直到到达一定次数或者满足自己所求的精度，停止迭代。
4、几何意义
第二步的公式中的f'(x0)可以看做△f(x0)/△x0=(f(x0)-0)/(x0-x1)=f(x0)/(x0-x1)，所以原式可化作x=x0-f(x0)/(f(x0/(x0-x1)))=x0-(x0-x1)=x1，这相当于做了一步什么操作呢，就是从(x0,f(x0))的点做函数f(x)的切线，交x轴于x1，然后从(x1,f(x1))的点做函数f(x)的切线，交x轴于x2，不断迭代，我们会发现f(xn)会越来越趋近于0，这也就是牛顿迭代的中心思想。
方法是有了，可是我们怎么用程序实现呢，因为公式中涉及了导数，需要用到微分，当然，你可以用数学的方式把导数的表达式都计算好了，然后再通过代码的形式带入，但是这显然不够智能，如果函数f(x)会发生变化的话，每次还要去改代码，非常的不灵活。
在此介绍一个非常好用的模块，python的sympy模块，通过这个模块，我们可以非常方便的算出微积分的结果表达式，下面来看demo：
| from sympy import *
# 定义一个数学符号x
x, y = symbols('x y')
# 定义一个数学表达式
exam = 3*x**4 + 5*x**3 + x**2 + 8*x + 6
# ======================求导
# 求导
print(exam.diff(x))
print(diff(exam, x))
# 结果：12*x**3 + 15*x**2 + 2*x + 8
# 求二次导
print(exam.diff(x, x))
print(diff(exam, x, x))
# 结果：2*(18*x**2 + 15*x + 1)
# 可以先将计算公式存储到一个变量中，然后再通过doit计算，其实和diff的效果是一样的
deriv = Derivative(exam, x, x)
print(deriv.doit())
# 结果：2*(18*x**2 + 15*x + 1)
# ======================求积分
# 求不定积分
print(integrate(exam, x))
print(exam.integrate(x))
# 结果：3*x**5/5 + 5*x**4/4 + x**3/3 + 4*x**2 + 6*x
# 二次不定积分
print(integrate(exam, x, x))
print(exam.integrate(x, x))
# 结果：x**6/10 + x**5/4 + x**4/12 + 4*x**3/3 + 3*x**2
# 定积分，x在0到2的积分
print(integrate(exam, (x, 0, 2)))
print(exam.integrate((x, 0, 2)))
# 结果：1048/15
# 二元积分（可定积分和不定积分结合）
print(integrate(x**2 + y**2, (x, 0, 1), (y, 0, 1)))
# 结果：2/3
print((x**2 + y**2).integrate((x, 0, 1), y))
# 结果：y**3/3 + y/3
# ======================求极限
print(limit(sin(x)/x, x, 0))
# 结果：1
print(limit(sin(x)/x, x, oo))
# 结果：0
# 当x趋于0时，从负方向或是正方向逼近，可能会有不同的结果
print(limit(1/x, x, 0, '+'))
# 结果：oo
print(limit(1/x, x, 0, '-'))
# 结果：-oo
# ======================求泰勒展开公式
# 传入的参数分别为x0的取值、泰勒展开项数
print((x**4).series(x, 2, 2))
# 结果：-48 + 32*x + O((x - 2)**2, (x, 2))
print(series((x**4), x, 2, 6))
# 结果：32*x + (x - 2)**4 + 8*(x - 2)**3 + 24*(x - 2)**2 - 48
|
好了，有了以上知识为基础就可以解决问题了，还是回到上面的公式，为了方便码字，我们假设小明每期的投资都是50元，上面的公式可表示为：
可以转化为f(x)=0的形式：
接下来初始化一个x0，比如0.2，然后一步步进行迭代x0，直到f(x0)足够接近0，或是达到一定轮数停止，即可得到近似解，下面用代码进行展示：
| from sympy import *
x = Symbol('x')
# 获取fx的表达式
def get_fx():
result = 10000
for i in range(1, 13):
result -= 50 * (1 + x)**i
return result
# 获取fx导数的表达式
def get_diff_fx():
fx = get_fx()
diff_fx = fx.diff(x)
return diff_fx
# 获取函数的值
def get_fx_value(fx, x, limit_val):
return limit(fx, x, limit_val)
# 牛顿迭代计算近似解
def newton_iter():
# 初值，根据经验选取，如果选取的值与最终求得值得接近，会减少迭代次数，提高程序效率
x0 = 0.2
# 精度要求，当误差小于此值时，结束迭代
prec = 0.0001
# 最多迭代20次
for i in range(15):
fx = get_fx_value(get_fx(), x, x0)
print(fx)
if abs(fx) <= prec:
break
diff_fx = get_fx_value(get_diff_fx(), x, x0)
x0 -= fx / diff_fx
print("共进行了%s次迭代，最终的解为%s，最终算得的fx的值为%s" % (i, x0, round(fx, 10)))
newton_iter()
# 运行结果：共进行了7次迭代，最终的解为0.403702825570710，最终算得的fx的值为-2.86490831058472e-11
|
以上为代码，我用x0=0.4又做了一次测试，发现只要迭代3次就是算出结果了，
| 共进行了3次迭代，最终的解为0.403702825570709，最终算得的fx的值为-2.86490831058472e-11
|
以上内容适用范围并不仅限于计算irr或是银行率，而是所有涉及到高阶方程的场景。
凡本网注明"来源：XXX "的文/图/视频等稿件，本网转载出于传递更多信息之目的，并不意味着赞同其观点或证实其内容的真实性。如涉及作品内容、版权和其它问题，请与本网联系，我们将在第一时间删除内容！
作者:
来源： https://www.jianshu.com/p/b4b0ccd8cf8c