Python - 实现数值方程求解器(Newton-Raphson) [英] Python - Implementing a numerical equation solver (Newton-Raphson)

问题描述

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

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

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

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

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

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.

1. 如何编写方程的代码?

2. 当近似值不再变化时，如何终止循环?

1. How do I write the code for the equation?

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

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

推荐答案

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

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...

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

输出:

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

编辑 - 您也可以检查 newY 的值并在它足够接近零"时停止 - 但通常您会保持这种状态直到 x 发生变化是 <=h(您可以争论数值方法中 = 符号的值 - 我自己更喜欢更强调的 <，认为再一次迭代不会有什么坏处.)

Edit - you could also check the value of newY and stop when it is "close enough to zero" - but usually you keep this going until the change in x is <=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.).

这篇关于Python - 实现数值方程求解器(Newton-Raphson)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！