制作牛顿法求函数根的动画

Question

在我的微积分课之前，我想创建一个图形，该图形会自我更新以显示牛顿法寻找函数根的过程。所以我的问题是，使用 Python 和 Sympy 来处理问题背后的所有混乱计算，我将如何以一种信息丰富的方式绘制这个图表。

乍一看，Sympy 的最新库似乎没有提供任何方法来简单地将多个图添加到现有图上，这是一个问题。

希望有人找到以下使用答案:)

Answer 1

然而，经过大约 2 天的工作，我设法想出了这个解决方案。我仍在使用 Sympy 作为导数计算和函数输入的后端，但我将这些表达式转换为 lambda 函数以生成可以输入 Matplotlib 坐标的值。我还清理了代码并将 cmets 添加到最重要的部分。

为了使用这个程序（仅在 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()