我警告你,这可能令人困惑,我写的代码更多的是思维导图而不是完成代码.
我正在尝试使用Newton-Raphson方法来求解方程.我无法弄清楚的是如何写这个
Python中的等式,用于计算最后一次近似(xn)的下一个近似值(xn + 1).我必须使用循环,越来越接近真实的答案,并且当近似值之间的变化小于变量h时,循环应该终止.
如何编写等式的代码?
当近似值不再变化时,如何终止循环?
def derivative(f, x, h): deriv = (1.0/(2*h))*(f(x+h)-f(x-h)) return deriv
def solve(f, x0, h): xn = x0 prev = 0 while ( approx - prev > h): xn = xn - (f(xn))/derivative(f, xn, h) return xn
Floris.. 8
以下是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
输出:
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
(你可以争论的价值=
符号的数值方法-我喜欢更加强调<
自己,认为再多一次迭代不会受到伤害.)
以下是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
输出:
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
(你可以争论的价值=
符号的数值方法-我喜欢更加强调<
自己,认为再多一次迭代不会受到伤害.)