制作牛顿法求函数根的动画
Making an animation of Newton's Method for finding roots of a function
在微积分 class 之前,我曾想过创建一个图表,该图表会自行更新以显示牛顿法求函数根的过程。所以我的问题是,使用 Python 和 Sympy 来处理问题背后的所有混乱演算,我将如何以信息丰富的方式绘制图表。
乍一看,Sympy 的最新库似乎没有提供任何方法来简单地将多个图添加到现有图上,这是一个问题。
希望有人找到以下有用的答案:)
经过大约 2 天的工作,我设法想出了这个解决方案。我仍在使用 Sympy 作为导数计算和函数输入的后端,但我正在将这些表达式转换为 lambda 函数以生成可以输入到 Matplotlib 坐标中的值。我还清理了代码并为最重要的部分添加了注释。
为了使用这个程序(只在 Windows 上测试过)我安装了最新的 Anaconda 包。
'''
Created on Aug 10, 2016
@author: Clement Hathaway
Example: enter in:
6*x**3 -73*x**2 -86*x +273
-20
-10
20
-2200
2000
GO!
'''
from numpy import linspace
from sympy.parsing.sympy_parser import parse_expr
from sympy import *
from sympy.solvers.solveset import solveset, solveset_real
import matplotlib.pyplot as plt
#from numpy.doc.constants import lines
def distanceFromRoot(x, values): ## Return the lowest distance between the supplied roots
currentDistance = abs(x-values[0])
if len(values) > 1: ## Only run if there is more than one root
for value in values:
newDist = abs(x-value)
if newDist < currentDistance:
currentDistance = newDist
return currentDistance
def generateTangents(x0, real_roots, accuracy):
if distanceFromRoot(x0, real_roots) > accuracy: ## Continue if the value we have calculated is still too far from a real root
tangent = simplify(dy.subs(x, x0)*(x-x0)+(y.subs(x,x0))) ## Find equation of tangent at our x0
tangent_intercept = solve(tangent, x)[0] ## Find tangent's x intercept
lam_tan = lambdify(x, tangent, modules=['numpy']) ## Function to generate y values of the tangent line
tan_y_vals = lam_tan(x_values) ## Y values of the tangent line to plot
text.set_text("Current x-estimate: " + str(N(tangent_intercept))) ## Display estimate
plt.plot(x_values, tan_y_vals, 'b-') ## Plot our new y values with a colour of blue and a ocnsistent line
plt.draw() ## Draw the plot
plt.pause(0.5) ## Pause for 1 second before drawing the next line
generateTangents(N(tangent_intercept), real_roots, accuracy) ## Call itself again recursively
else: ## Do nothing
## Display some ending text
print("Found value to be: " + str(x0))
print("With real roots: " + str(real_roots))
pass
if __name__ == '__main__':
## Equation
x = symbols('x')
y = parse_expr(str(raw_input("Enter Equation: ")))
dy = diff(y, x) ## Differentiate y in terms of x
x0 = float(input("Enter starting value for x: "))
roots = solveset_real(y, x) ## Find roots of y in terms of x
roots_array = [] ## Get in terms of array
for root in roots:
roots_array.append(N(root))
## Setup Values
x_min = int(input("Enter x-axis min: "))
x_max = int(input("Enter x-axis max: "))
y_min = int(input("Enter y-axis min: "))
y_max = int(input("Enter y-axis max: "))
res = 200
## Convert to use with other library
lam_y = lambdify(x, y, modules=['numpy']) ## Functions of our main func
lam_dy = lambdify(x, dy, modules=['numpy']) ## Differential of our eq
## Calculate starting values
x_values = linspace(x_min, x_max, res) ## Array of x values between xmin and max seperated by res
y_values = lam_y(x_values) ## Calculated y values
## Graph Setup
plt.axis([x_min, x_max, y_min, y_max]) ## Setup graph axis
plt.grid() ## Enable the graph grid
plt.ion() ## Make graph interactive to add future lines
## Plot main Function
y0_values = linspace(0,0, res)
plt.plot(x_values, y_values, 'g-', linewidth=2)
plt.plot(x_values, y0_values,'k-', linewidth=2) ## Create more noticable x axis
text = plt.text(0,5,"Current x-estimate: " + str(x0), fontsize = 30)
plt.draw()
plt.pause(0.1)
## Start generating tangent lines! weo
generateTangents(x0, roots_array, 0.0000001)
## Keep the graph from disappearing
plt.ioff()
plt.show()
在微积分 class 之前,我曾想过创建一个图表,该图表会自行更新以显示牛顿法求函数根的过程。所以我的问题是,使用 Python 和 Sympy 来处理问题背后的所有混乱演算,我将如何以信息丰富的方式绘制图表。
乍一看,Sympy 的最新库似乎没有提供任何方法来简单地将多个图添加到现有图上,这是一个问题。
希望有人找到以下有用的答案:)
经过大约 2 天的工作,我设法想出了这个解决方案。我仍在使用 Sympy 作为导数计算和函数输入的后端,但我正在将这些表达式转换为 lambda 函数以生成可以输入到 Matplotlib 坐标中的值。我还清理了代码并为最重要的部分添加了注释。
为了使用这个程序(只在 Windows 上测试过)我安装了最新的 Anaconda 包。
'''
Created on Aug 10, 2016
@author: Clement Hathaway
Example: enter in:
6*x**3 -73*x**2 -86*x +273
-20
-10
20
-2200
2000
GO!
'''
from numpy import linspace
from sympy.parsing.sympy_parser import parse_expr
from sympy import *
from sympy.solvers.solveset import solveset, solveset_real
import matplotlib.pyplot as plt
#from numpy.doc.constants import lines
def distanceFromRoot(x, values): ## Return the lowest distance between the supplied roots
currentDistance = abs(x-values[0])
if len(values) > 1: ## Only run if there is more than one root
for value in values:
newDist = abs(x-value)
if newDist < currentDistance:
currentDistance = newDist
return currentDistance
def generateTangents(x0, real_roots, accuracy):
if distanceFromRoot(x0, real_roots) > accuracy: ## Continue if the value we have calculated is still too far from a real root
tangent = simplify(dy.subs(x, x0)*(x-x0)+(y.subs(x,x0))) ## Find equation of tangent at our x0
tangent_intercept = solve(tangent, x)[0] ## Find tangent's x intercept
lam_tan = lambdify(x, tangent, modules=['numpy']) ## Function to generate y values of the tangent line
tan_y_vals = lam_tan(x_values) ## Y values of the tangent line to plot
text.set_text("Current x-estimate: " + str(N(tangent_intercept))) ## Display estimate
plt.plot(x_values, tan_y_vals, 'b-') ## Plot our new y values with a colour of blue and a ocnsistent line
plt.draw() ## Draw the plot
plt.pause(0.5) ## Pause for 1 second before drawing the next line
generateTangents(N(tangent_intercept), real_roots, accuracy) ## Call itself again recursively
else: ## Do nothing
## Display some ending text
print("Found value to be: " + str(x0))
print("With real roots: " + str(real_roots))
pass
if __name__ == '__main__':
## Equation
x = symbols('x')
y = parse_expr(str(raw_input("Enter Equation: ")))
dy = diff(y, x) ## Differentiate y in terms of x
x0 = float(input("Enter starting value for x: "))
roots = solveset_real(y, x) ## Find roots of y in terms of x
roots_array = [] ## Get in terms of array
for root in roots:
roots_array.append(N(root))
## Setup Values
x_min = int(input("Enter x-axis min: "))
x_max = int(input("Enter x-axis max: "))
y_min = int(input("Enter y-axis min: "))
y_max = int(input("Enter y-axis max: "))
res = 200
## Convert to use with other library
lam_y = lambdify(x, y, modules=['numpy']) ## Functions of our main func
lam_dy = lambdify(x, dy, modules=['numpy']) ## Differential of our eq
## Calculate starting values
x_values = linspace(x_min, x_max, res) ## Array of x values between xmin and max seperated by res
y_values = lam_y(x_values) ## Calculated y values
## Graph Setup
plt.axis([x_min, x_max, y_min, y_max]) ## Setup graph axis
plt.grid() ## Enable the graph grid
plt.ion() ## Make graph interactive to add future lines
## Plot main Function
y0_values = linspace(0,0, res)
plt.plot(x_values, y_values, 'g-', linewidth=2)
plt.plot(x_values, y0_values,'k-', linewidth=2) ## Create more noticable x axis
text = plt.text(0,5,"Current x-estimate: " + str(x0), fontsize = 30)
plt.draw()
plt.pause(0.1)
## Start generating tangent lines! weo
generateTangents(x0, roots_array, 0.0000001)
## Keep the graph from disappearing
plt.ioff()
plt.show()