使用牛顿迭代法求高次方程的根

使用牛顿迭代法求高次方程的根

引用

CSDN

1.

https://blog.csdn.net/weixin_59908804/article/details/143134981

牛顿迭代法是一种用于求解方程近似根的数值方法，其基本思想是通过不断求曲线上的某点的切线，将切线的零点近似为方程的解，如此循环迭代。本文将详细介绍牛顿迭代法的原理，并通过一个金融领域的实际案例（计算内部收益率IRR）来展示其具体应用。

牛顿迭代法

牛顿迭代法通过不断求曲线上的某点的切线，将切线的零点近似为方程的解，如此循环迭代。牛顿迭代法以平方速度逼近函数的零点，是比二分法更快的迭代方法。

数学解释

对于方程f(x) = 0，取定一个r0，求f(x)在r0处的切线l。然后求l的零点r1，然后以此迭代，依次求r2, r3, ...

实例展示

在金融中，我们有时会用内部收益率IRR来评价项目的投资财务效益，它等于使得投资净现值NPV等于0的贴现率。换句话说，给定项目的期数T、初始现金流CF0和项目各期的现金流CF1, CF2, …，CFT，IRR是下面方程的解：

为了简单起见，本题假定：除了项目启动时有一笔投入（即初始现金流CF0 < 0）之外，其余各期均能赚钱（即对于所有i=1,2,…,T，CFi > 0）。根据定义，IRR可以是负数，但始终大于-1。

题目分析

本题本质上是求方程NPV(IRR+1) = 0的解。注意到，实际上IRR+1总是作为一个整体，所以我们记x = IRR+1。故而求NPV(x)（x>0）的零点即可。

不难发现，NPV'(x)始终小于0，所以NPV(x)最多只有一个零点。

求ri处的切线：
y - NPV(x) = NPV'(x)(x - ri)
迭代：ri+1 = -NPV(ri) / NPV'(ri) + ri

此外，对于r0的给出，根据经验给出1.5（具体是怎么知道的我也不清楚）。由此反复迭代，直到NPV(x)符合条件，此时以rn近似x0。

代码实现

求NPV(x)：

double functionF(double x,double*CF,int len) {
    double sum = 0.0;
    for (int i = 0; i < len; i++)
        sum += CF[i] * pow(x, -i);
    return sum;
}

求NPV'(x)：

double derivitiveF(double x, double* CF, int len) {
    double sum = 0.0;
    for (int i = 1; i < len; i++)
        sum += -CF[i] * pow(x, -i - 1);
    return sum;
}

迭代：

NPV = functionF(x, CF, T + 1);
if (NPV < 1e-5)break;//退出条件
x = x - NPV / derivitiveF(x, CF, T + 1);

完整代码：

#define MAX_COUNT 20
double functionF(double x,double*CF,int len) {
    double sum = 0.0;
    for (int i = 0; i < len; i++)
        sum += CF[i] * pow(x, -i);
    return sum;
}
double derivitiveF(double x, double* CF, int len) {
    double sum = 0.0;
    for (int i = 1; i < len; i++)
        sum += -CF[i] * pow(x, -i - 1);
    return sum;
}
int main(){
    float result[MAX_COUNT];
    int index = 0;
    while(1){//反复接受数据 知道接收到T=0退出循环
        unsigned int T = 0;
        double CF[MAX_COUNT];//storage the CFn
        double NPV = 0.0;
        double x = 1.5;//initial value
        scanf_s("%d", &T);
        if (0 == T)break;
        printf("You should print at least %d datas.\n", T + 1);
        for (int i = 0; i < T + 1; i++)scanf_s("%lf", &CF[i]);
        while (1) {
            NPV = functionF(x, CF, T + 1);
            if (NPV < 1e-5)break;
            x = x - NPV / derivitiveF(x, CF, T + 1);//interate
        }
        double IRF = x - 1;
        result[index++] = (float)IRF;//storage
    }
    for(int i=0;i<index;i++)//print the results
        printf("%0.2f\n",result[i]);
    return 0;
}