Python - 实现数值方程求解器(Newton-Raphson)

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/20659456/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me): StackOverFlow

提示:将鼠标放在中文语句上可以显示对应的英文。显示中英文
时间:2020-08-18 20:58:36  来源:igfitidea点击:

Python - Implementing a numerical equation solver (Newton-Raphson)

pythonwhile-loopnumerical-methods

提问by user2906011

I am warning you, this could be confusing, and the code i have written is more of a mindmap than finished code..

我警告你,这可能会令人困惑,我写的代码更像是一个思维导图而不是完成的代码。

I am trying to implement the Newton-Raphson method to solve equations. What I can't figure out is how to write this

我正在尝试实施 Newton-Raphson 方法来求解方程。我想不通的是如何写这个

enter image description here

在此处输入图片说明

equation in Python, to calculate the next approximation (xn+1) from the last approximation (xn). I have to use a loop, to get closer and closer to the real answer, and the loop should terminate when the change between approximations is less than the variable h.

Python 中的方程,从上一个近似值 (xn) 计算下一个近似值 (xn+1)。我必须使用循环来越来越接近真正的答案,并且当近似值之间的变化小于变量 h 时,循环应该终止。

  1. How do I write the code for the equation?
  2. How do I terminate the loop when the approximations are not changing anymore?

    Calculates the derivative for equation f, in point x, with the accuracy h (this is used in the equation for solve())

    def derivative(f, x, h):
        deriv = (1.0/(2*h))*(f(x+h)-f(x-h))
        return deriv
    

    The numerical equation solver

    Supposed to loop until the difference between approximations is less than h

    def solve(f, x0, h):
        xn = x0
        prev = 0
    
        while ( approx - prev > h):
             xn = xn - (f(xn))/derivative(f, xn, h)
    
        return xn
    
  1. 我如何编写方程的代码?
  2. 当近似值不再改变时,如何终止循环?

    计算方程 f 的导数,在点 x,精度为 h(这用于求解()的方程)

    def derivative(f, x, h):
        deriv = (1.0/(2*h))*(f(x+h)-f(x-h))
        return deriv
    

    数值方程求解器

    应该循环直到近似值之间的差异小于 h

    def solve(f, x0, h):
        xn = x0
        prev = 0
    
        while ( approx - prev > h):
             xn = xn - (f(xn))/derivative(f, xn, h)
    
        return xn
    

回答by Bathsheba

You'll be lucky to get convergence since your derivative is not exact to the limits of numerical precision.

你会很幸运地得到收敛,因为你的导数不精确到数值精度的限制。

Aside from your not guarding against division by zero, there's nothing wrong with your implementation of the Newton-Raphson result: that is your statement xn = xn - (f(xn))/derivative(f, xn, h)

除了您没有防止被零除之外,您对 Newton-Raphson 结果的实现没有任何问题:这就是您的陈述 xn = xn - (f(xn))/derivative(f, xn, h)

But since you're using an approximate derivative, you should switch to a different optimisation scheme once you have the root bracketed. So as far as the Newton Raphson part is concerned, thatis your termination condition. A good optimiser to use is brentwhich will always find a root once bracketed; even without a derivative.

但是,由于您使用的是近似导数,因此在将根括号括起来后,您应该切换到不同的优化方案。因此,就 Newton Raphson 部分而言,就是您的终止条件。一个很好的优化器是brent,它总是会在括号中找到一个根;即使没有衍生品。

回答by Ramchandra Apte

How do I write the code for the equation?

我如何编写方程的代码?

The Newton-Raphson method actuallyfinds the zeroes of a function. To solve an equation g(x) = y, one has to make the function passed to the solver g(x)-yso that when the function passed to the solver gives zero, g(x)=y.

Newton-Raphson 方法实际上是找到函数的零点。要求解方程g(x) = y,必须将函数传递给求解器,g(x)-y以便当传递给求解器的函数为零时,g(x)=y

How do I terminate the loop when the approximations are not changing anymore?

当近似值不再改变时,如何终止循环?

Your code already does that as when the two approximations are equal to each other, the difference will be 0, which is less than h.

您的代码已经这样做了,因为当两个近似值彼此相等时,差异将为0,小于h

current_approx - previous_approx > hshould be current_approx - previous_approx >= hsince you want it to terminate when the difference is less than h. Also improved the variable names.

current_approx - previous_approx > h应该是current_approx - previous_approx >= h因为您希望它在差异小于 时终止h。还改进了变量名称。

def derivative(f, x, accuracy):
    deriv = (1.0/(2*accuracy))*(f(x+accuracy)-f(x-accuracy))
    return deriv

def solve(f, approximation, h):
    current_approx = approximation
    prev_approximation = float("inf") # So that the first time the while loop will run

    while current_approx - prev_approx >= h:
        prev_approximation = current_approx
        current_approx = current_approx - (f(current_approx))/derivative(f, current_approx, h)

    return current_approx

回答by Floris

Here is the implementation of a N-R solver expanding what you wrote above (complete, working). I added a few extra lines to show what is happening...

这是扩展你上面写的内容的 NR 求解器的实现(完成,工作)。我添加了一些额外的行来显示正在发生的事情......

def derivative(f, x, h):
      return (f(x+h) - f(x-h)) / (2.0*h)  # might want to return a small non-zero if ==0

def quadratic(x):
    return 2*x*x-5*x+1     # just a function to show it works

def solve(f, x0, h):
    lastX = x0
    nextX = lastX + 10* h  # "different than lastX so loop starts OK
    while (abs(lastX - nextX) > h):  # this is how you terminate the loop - note use of abs()
        newY = f(nextX)                     # just for debug... see what happens
        print "f(", nextX, ") = ", newY     # print out progress... again just debug
        lastX = nextX
        nextX = lastX - newY / derivative(f, lastX, h)  # update estimate using N-R
    return nextX

xFound = solve(quadratic, 5, 0.01)    # call the solver
print "solution: x = ", xFound        # print the result

output:

输出:

f( 5.1 ) =  27.52
f( 3.31298701299 ) =  6.38683083151
f( 2.53900845771 ) =  1.19808560807
f( 2.30664271935 ) =  0.107987672721
f( 2.28109300639 ) =  0.00130557566462
solution: x =  2.28077645501

Edit - you could also check the value of newYand stop when it is "close enough to zero" - but usually you keep this going until the change in xis <=h(you can argue about the value of the =sign in a numerical method - I prefer the more emphatic <myself, figuring that one more iteration won't hurt.).

编辑 - 你也可以检查 的值newY并在它“足够接近零”时停止 - 但通常你会一直这样做直到改变x<=h(你可以争论=数字方法中符号的值- 我更喜欢<我自己更加强调,认为再一次迭代不会受到伤害。)。

回答by 74U n3U7r1no

If code under 'try:' cannot be implemented and the compiler is given the error 'ZeroDivisionError' then it will execute the code under 'except ZeroDivisionError:'. Although, if you'd like to account for another compiler exception 'XYZ' with a specific code implementation then add an additional 'except XYZ:"

如果“try:”下的代码无法实现并且编译器给出错误“ZeroDivisionError”,那么它将执行“除了 ZeroDivisionError:”下的代码。虽然,如果您想用特定的代码实现来解释另一个编译器异常“XYZ”,然后添加一个额外的“除了 XYZ”:

 try:
    nextX = lastX - newY / derivative(function, lastX, h)
 except ZeroDivisionError:
    return newY