牛顿法求解方程（python和C++）

系统管理员 2022-05-21 00:19 300阅读 0赞

# 牛顿法 #

> **GitHub: [https://github.com/Sean16SYSU/Algorithms4N][https_github.com_Sean16SYSU_Algorithms4N]**

目标函数为:  f ( x ) = 0 f(x) = 0 f(x)=0

算法很简单如下:  
 x n + 1 = x n − f ( x n ) f ′ ( x n ) x\_\{n+1\} = x\_\{n\} - \\frac\{f(x\_n)\}\{f'(x\_n)\} xn\+1=xn−f′(xn)f(xn)

其实是对f(x)做了泰勒展开, 用前两项来做对于f(x)的近似替代  
 f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) 2 ! ( x − x 0 ) 2 + . . . f(x) = f(x\_0) + f'(x\_0)(x-x\_0) + \\frac\{f''(x\_0)\}\{2!\}(x-x\_0)^2+... f(x)=f(x0)\+f′(x0)(x−x0)\+2!f′′(x0)(x−x0)2\+...

即,  f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) = 0 f(x\_0) + f'(x\_0)(x-x\_0) = 0 f(x0)\+f′(x0)(x−x0)=0

故有,  x n + 1 = x n − f ( x n ) f ′ ( x n ) x\_\{n+1\} = x\_\{n\} - \\frac\{f(x\_n)\}\{f'(x\_n)\} xn\+1=xn−f′(xn)f(xn)

# Newton-Raphson method #

Objective function:  f ( x ) = 0 f(x) = 0 f(x)=0

The simple but useful iteration formula:  
 x n + 1 = x n − f ( x n ) f ′ ( x n ) x\_\{n+1\} = x\_\{n\} - \\frac\{f(x\_n)\}\{f'(x\_n)\} xn\+1=xn−f′(xn)f(xn)

In fact, the first two terms of the Taylor Expansion of f(x) are ultilized to make an approximately subsitution of f(x)  
 f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) 2 ! ( x − x 0 ) 2 + . . . f(x) = f(x\_0) + f'(x\_0)(x-x\_0) + \\frac\{f''(x\_0)\}\{2!\}(x-x\_0)^2+... f(x)=f(x0)\+f′(x0)(x−x0)\+2!f′′(x0)(x−x0)2\+...

namely,  f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) = 0 f(x\_0) + f'(x\_0)(x-x\_0) = 0 f(x0)\+f′(x0)(x−x0)=0

So,  x n + 1 = x n − f ( x n ) f ′ ( x n ) x\_\{n+1\} = x\_\{n\} - \\frac\{f(x\_n)\}\{f'(x\_n)\} xn\+1=xn−f′(xn)f(xn)

上述的迭代函数，本质上，其实是关于 f ( x ) f(x) f(x)的微分中值定理的变形，也可以说是对于一阶的泰勒展开的变形。其实，是一种切线下降法，对于低精度下的要求的话，下降应该是非常迅速的。

## 代码 ##

*  c++实现根号计算

#include <iostream>
    using namespace std;
    #define itTime 1000 // the Maximum Iteration Times 
    #define eps 1e-10 // precision accuracy
    #define abs(x) ((x > 0? (x):-(x))) 
    int main() { 
    	int N, time = 1;
    	cin >> N; // input number N;
    	double x = N / 2., tmp;
    	while (time++ < itTime) { 
    		tmp = x / 2. + N / (2. * x);
    		if (abs(tmp - x) < eps) break;
    		x = tmp;
    	}
    	printf("sqrt(%d) = %.10lf\n", N, tmp);
    	system("pause");
    }

*  Python计算

from sympy import *
    import random
    
    x = symbols("x")
    
    func = x*exp(x)-1
    ffunc = diff(func, x)
    
    begin = 1
    end = 2
    
    MAXSTEP = 100
    
    step_count = 0
    
    x0 = random.uniform(begin, end)
    temp = func.subs(x, x0)
    
    while step_count < MAXSTEP and abs(temp) > 1e-10:
        x0 = x0 - (temp / (ffunc.subs(x, x0)))
        temp = func.subs(x, x0)
        step_count += 1
    print(x0)
    print(step_count)

[https_github.com_Sean16SYSU_Algorithms4N]: https://github.com/Sean16SYSU/Algorithms4N