用python算微积分及牛顿迭代求解高阶方程

前段时间,工作中遇到了一个难题,网上查了很多资料也没有明确的解决方案,后来自己想到一种解决方法,决定分享出来。
下面先说下问题:如何计算irr(内部收益率)
我们先来看下百度百科对这1名词的解释和计算方式的描述:


image.png

image.png

看完上面的解释,给了一大堆公式,还是没有告诉怎么求,查看其它资料也大部分言尽于此,固执的我看了更多的名词解释,终于明白了这个是怎么来的,大概就是下面这个公式:


image.png

上面公式中的y和a1到an都是已知量,x就是要求的irr。这个公式怎么来的,或者说在生活中的具体例子是什么呢,可以这样理解:
小明往银行存钱,假设月利率是1%(只为举例,实际生活中肯定是低于此值的),1月到12月每月存入100块,问到年底小明能拿回多少钱,肯定是像下图这样算。


image.png

那么问题再升级一下,如果小明每月存入的金额不固定了,用an表示每月存入的金额,那结果怎么算呢?那其实还是很简单,把上面公式的100改成a1到a12就好了,如下公式:


image.png

那么问题再次升级,如果告诉你小明年底能拿到1000元,每期存的金额分别是a1到a12,假设利率是稳定不变的,此时问你利率是多少,你该怎么算呢?首先肯定还是列出下面的公式:


image.png

公式是有了,但是这是一个高阶方程,该如何求解呢,纵观数学史上,方程一旦高于5次是没有求根公式的,但是伟大的先贤们还是找到了一些解决方案,就是通过多次迭代之后,近似的找到方程的解。今天要介绍的就是几种迭代方案中的一种——牛顿迭代,用于求解f(x)=0的解:

下面我们来看看牛顿迭代的原理:
首先我们根据泰勒公式,我们知道f(x)可由其邻域x0的n阶导数的表达式近似表示,公式如下:


image.png

因为公式中后面项的分母越来越大,其对f(x)值的影响就越来越小,所以我们取线性部分,再通过迭代的方式可达到近似求解的目的:
1、取泰勒展开的线性部分


image.png

2、解线性部分的方程


image.png

image.png

3、对第2部不断进行迭代,直到到达一定次数或者满足自己所求的精度,停止迭代。
4、几何意义


image.png

第二步的公式中的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元,上面的公式可表示为:


image.png

可以转化为f(x)=0的形式:


image.png

接下来初始化一个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或是银行率,而是所有涉及到高阶方程的场景。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 202,905评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,140评论 2 379
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,791评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,483评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,476评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,516评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,905评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,560评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,778评论 1 296
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,557评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,635评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,338评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,925评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,898评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,142评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,818评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,347评论 2 342

推荐阅读更多精彩内容