使用 sympy 逐步微分

2024-04-07

我正在尝试制作一个 python proram 来查找导数和积分以及展示如何进行。到目前为止,我发现有一个integral_steps函数返回所使用的步骤,但我还没有找到等效的微分函数。

有谁知道是否有等效的?

如果没有,您对如何找到找到导数所需的步骤有什么想法吗?


方法一(手动)

查看代码,派生类 https://github.com/sympy/sympy/blob/723630c/sympy/core/function.py#L1048是顶层逻辑所在的地方。这只是顶层部分。从那时起,计算需要计算内部不同节点的导数表达式树 https://en.wikipedia.org/wiki/Abstract_syntax_tree.

表达式树的每个特定节点的逻辑位于_eval_derivative method https://github.com/sympy/sympy/search?q=%22_eval_derivative%22对应于每个特定的节点类型。

这将允许您向这些添加代码_eval_derivative方法以便追踪整个过程并找到所有步骤。

方法2(使用示踪剂)

Python有多个tracing https://en.wikipedia.org/wiki/Tracing_(software)包。蟒蛇猎人 https://python-hunter.readthedocs.org写于@ionelmc https://stackoverflow.com/users/23658/ionelmc实际上非常好并且非常适合这个用例。

在许多其他功能中,它允许在函数开始执行时安装某些回调,并在函数返回其值时安装另一个回调。事实上这正是我们所需要的。

下面是一个展示如何使用它的示例(我在 Python 3.7.3、SymPy 1.7 和 Hunter 3.3.1 上运行并测试了它):

import hunter
import sys 
from hunter import Q, When, Stop
hunter.trace(
        Q(module_contains="sympy",function='_eval_derivative',kind_in=["call","return"],action=hunter.CallPrinter(repr_func=str))
        )


from sympy import *
x = symbols('x')
f = 1/(x * sin(x)**2)
f.diff(x)

因此,这允许我们选择我们想要检查的数据结构,我们想要如何打印它们,并且它允许我们看到微分过程的中间步骤:

[...]7/site-packages/sympy/core/power.py:1267  call      => _eval_derivative(self=sin(x)**(-2), s=x)
[...]7/site-packages/sympy/core/power.py:1267  call         => _eval_derivative(self=<sympy.core.power.Pow object at 0x7f5925337150>, s=<sympy.core.symbol.Symbol object at 0x7f5925b6a2b0>)
[...]ite-packages/sympy/core/function.py:598   call            => _eval_derivative(self=sin(x), s=x)
[...]ite-packages/sympy/core/function.py:598   call               => _eval_derivative(self=<sympy.functions.elementary.trigonometric.sin object at 0x7f592589ee08>, s=<sympy.core.symbol.Symbol object at 0x7f5925b6a2b0>)
[...]ite-packages/sympy/core/function.py:612   return             <= _eval_derivative: cos(x)
[...]ite-packages/sympy/core/function.py:612   return          <= _eval_derivative: <sympy.functions.elementary.trigonometric.cos object at 0x7f592525fef8>
[...]7/site-packages/sympy/core/power.py:1271  return       <= _eval_derivative: -2*cos(x)/sin(x)**3
[...]7/site-packages/sympy/core/power.py:1271  return    <= _eval_derivative: <sympy.core.mul.Mul object at 0x7f5925259b48>
[...]7/site-packages/sympy/core/power.py:1267  call      => _eval_derivative(self=1/x, s=x)
[...]7/site-packages/sympy/core/power.py:1267  call         => _eval_derivative(self=<sympy.core.power.Pow object at 0x7f5925337200>, s=<sympy.core.symbol.Symbol object at 0x7f5925b6a2b0>)
[...]7/site-packages/sympy/core/power.py:1271  return       <= _eval_derivative: -1/x**2
[...]7/site-packages/sympy/core/power.py:1271  return    <= _eval_derivative: <sympy.core.mul.Mul object at 0x7f5925259f10>

如果您还想覆盖diff函数你可以改变上面的代码并有function_in=['_eval_derivative','diff']。这样不仅可以看到部分结果,还可以看到调用的部分diff函数及其返回值。

方法 3(使用跟踪器,构建调用图并将其可视化)

使用 graphviz、latex 和示踪剂(同样,蟒蛇猎人 https://python-hunter.readthedocs.io/)你实际上可以看到调用图 https://en.wikipedia.org/wiki/Call_graph更清楚。渲染每个中间步骤的所有公式确实需要一些时间,因为pdflatex正在使用(不过我确信乳胶有更快的渲染器)。

每个节点的值采用以下格式:

function_name
argument => return_value

好像有几个diff参数等于返回值的节点,我目前不确定如何解释。

如果该图以某种方式提到了每个规则的应用位置(我想不出一个简单的方法来做到这一点),那么它可能会更有用。

这也是代码:

import hunter
import sys
from hunter import Q, When, Stop, Action
from hunter.actions import  ColorStreamAction

formula_ltx = r'''
\documentclass[border=2pt,varwidth]{letter}
\usepackage{amsmath}
\pagenumbering{gobble}
\begin{document}
\[ \texttt{TITLE} \]
\[ FORMULA \]
\end{document}
'''

# ==============
# == Tracing ===
# ==============

from sympy.printing.latex import LatexPrinter, print_latex, latex

global call_tree_root

# a node object to hold an observed function call
# with its argument, its return value and its function name
class Node(object):
    def __init__(self, arg=None, retval=None, func_name=None):
        self.arg = arg
        self.retval = retval
        self.arg_ascii = ""
        self.retval_ascii = ""
        self.func_name = func_name
        self.uid = 0
        self.children = []

# this is a hunter action where we build a call graph and populate it
# so we can later render it
#
# CGBAction (Call Graph Builder Action)
class CGBAction(ColorStreamAction):
    def __init__(self, *args, **kwargs):
        super(ColorStreamAction, self).__init__(*args, **kwargs)
        # a custom call stack
        self.tstack = []
        global call_tree_root
        call_tree_root = Node(arg="",func_name="root")
        self.node_idx = 1
        self.tstack.append(call_tree_root)

    def __call__(self, event):
        if event.kind in ['return','call']:
            if event.kind == 'return':
                print(str(event.arg))
                if len(self.tstack) > 0:
                    top = self.tstack.pop()
                    top.retval = latex(event.arg)
                    top.retval_ascii = str(event.arg)

            elif event.kind == 'call':
                print(str(event.locals.get('self')))
                new = Node()
                new.uid = self.node_idx
                new.arg = latex(event.locals.get('self'))
                new.arg_ascii = str(event.locals.get('self'))
                top = self.tstack[-1]
                self.tstack.append(new)
                top.children.append(new)
                new.func_name = event.module + ":" + event.function
                self.node_idx += 1

hunter.trace(
        Q(module_contains="sympy",function_in=['_eval_derivative','diff'],kind_in=["call","return"],action=CGBAction)
        )

from sympy import *
x = symbols('x')
f = 1 / (x * sin(x)**2)
#f = 1 / (x * 3)
#f = sin(exp(cos(x)*asin(x)))
f.diff(x)

# ============================
# == Call graph rendering ====
# ============================

import os
import re
OUT_DIR="formulas"

if not os.path.exists(OUT_DIR):
    os.mkdir(OUT_DIR)

def write_formula(prefix,uid,formula,title):
    TEX = uid + prefix + ".tex"
    PDF = uid + prefix + ".pdf"
    PNG = uid + prefix + ".png"

    TEX_PATH = OUT_DIR + "/" + TEX
    with open(TEX_PATH,"w") as f:
        ll = formula_ltx
        ll = ll.replace("FORMULA",formula)
        ll = ll.replace("TITLE",title)
        f.write(ll)

    # compile formula
    CMD = """
        cd formulas ; 
        pdflatex {TEX} ;
        convert -trim -density 300 {PDF} -quality 90 -colorspace RGB {PNG} ;
    """.format(TEX=TEX,PDF=PDF,PNG=PNG)

    os.system(CMD)

buf_nodes = ""
buf_edges = ""
def dfs_tree(x):
    global buf_nodes, buf_edges

    arg = ("" if x.arg is None else x.arg)
    rv  = ("" if x.retval is None else x.retval)
    arg = arg.replace("\r","")
    rv = rv.replace("\r","")

    formula = arg + "\\Rightarrow " + rv
    print(x.func_name + " -> " + x.arg_ascii + " -> " + x.retval_ascii)

    x.func_name = x.func_name.replace("_","\\_")
    write_formula("",str(x.uid),formula,x.func_name)

    buf_nodes += """
    {0} [image="{0}.png" label=""];
    """.format(x.uid);

    for y in x.children:
        buf_edges += "{0} -> {1};\n".format(x.uid,y.uid);
        dfs_tree(y)

dfs_tree(call_tree_root)

g = open(OUT_DIR + "/graph.dot", "w")
g.write("digraph g{")
g.write(buf_nodes)
g.write(buf_edges)
g.write("}\n")
g.close()
os.system("""cd formulas ; dot -Tpng graph.dot > graph.png ;""")

将 SymPy 逻辑映射到微分规则

我认为剩下的一步是将中间节点从 SymPy 映射到微分规则 https://en.wikipedia.org/wiki/Differentiation_rules。以下是我能够绘制的一些内容:

  • 产品规则 https://en.wikipedia.org/wiki/Product_rule映射到sympy.core.mul.Mul._eval_derivative https://github.com/sympy/sympy/blob/3c45367/sympy/core/mul.py#L905
  • 链式法则 https://en.wikipedia.org/wiki/Chain_rule映射到sympy.core.function.Function._eval_derivative https://github.com/sympy/sympy/blob/3c453674cf39d18c70c27bb20d074ccd248304fe/sympy/core/function.py#L598
  • Sum rule https://en.wikipedia.org/wiki/Differentiation_rules#Differentiation_is_linear映射到sympy.core.add.Add._eval_derivative https://github.com/sympy/sympy/blob/3c453674cf39d18c70c27bb20d074ccd248304fe/sympy/core/add.py#L403
  • 求和的导数映射为sympy.concrete.summations.Sum._eval_derivative https://github.com/sympy/sympy/blob/706007ca2fe279020e099d36dd1db0e33123ac4c/sympy/concrete/summations.py#L267
  • 乘积的 n 阶导数莱布尼兹一般法则 https://en.wikipedia.org/wiki/General_Leibniz_rule映射到sympy.core.mul.Mul._eval_derivative_n_times https://github.com/sympy/sympy/blob/26703f859341dc06f00be81247083634671c4220/sympy/core/mul.py#L985
  • 广义幂律 https://en.wikipedia.org/wiki/Differentiation_rules#Generalized_power_rule映射到sympy.core.power.Pow._eval_derivative https://github.com/sympy/sympy/blob/706007ca2fe279020e099d36dd1db0e33123ac4c/sympy/core/power.py#L1269

我还没见过分数课sympy.core所以也许商规则 https://en.wikipedia.org/wiki/Quotient_rule是通过间接处理产品规则 https://en.wikipedia.org/wiki/Product_rule, and a 广义幂律 https://en.wikipedia.org/wiki/Differentiation_rules#Generalized_power_rule指数为-1。

Running

为了让它运行,你需要:

sudo apt-get install graphviz imagemagick texlive texlive-latex-base

还有文件/etc/ImageMagick-6/policy.xml必须使用以下行进行更新,以允许从PDF->PNG:

<policy domain="coder" rights="read|write" pattern="PDF" />

还有另一个调用图库叫做jonga https://github.com/bwohlberg/jonga但它有点通用,并且不允许完全过滤掉不需要的呼叫。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

使用 sympy 逐步微分 的相关文章

  • InterfaceError:连接已关闭(使用 django + celery + Scrapy)

    当我在 Celery 任务中使用 Scrapy 解析函数 有时可能需要 10 分钟 时 我得到了这个信息 我用 姜戈 1 6 5 django celery 3 1 16 芹菜 3 1 16 psycopg2 2 5 5 我也使用了psyc
  • 如何限制 sympy FiniteSet 包含符号

    我对 sympy 还很陌生 我尝试使用 linsolve 求解线性方程组 这产生了一个可以用以下两行重现的解决方案 d symbols d solution sets FiniteSet d 1 d 4 d 5 d 我的解决方案遵循限制 即
  • DreamPie 不适用于 Python 3.2

    我最喜欢的 Python shell 是DreamPie http dreampie sourceforge net 我想将它与 Python 3 2 一起使用 我使用了 添加解释器 DreamPie 应用程序并添加了 Python 3 2
  • 导入错误:没有名为 _ssl 的模块

    带 Python 2 7 的 Ubuntu Maverick 我不知道如何解决以下导入错误 gt gt gt import ssl Traceback most recent call last File
  • 打破嵌套循环[重复]

    这个问题在这里已经有答案了 有没有比抛出异常更简单的方法来打破嵌套循环 在Perl https en wikipedia org wiki Perl 您可以为每个循环指定标签 并且至少继续一个外循环 for x in range 10 fo
  • 从列表中的数据框列中搜索部分字符串匹配 - Pandas - Python

    我有一个清单 things A1 B2 C3 我有一个 pandas 数据框 其中有一列包含用分号分隔的值 某些行将包含与上面列表中的一项的匹配 它不会是完美的匹配 因为它在其中包含字符串的其他部分 该列 例如 该列中的一行可能有 哇 这里
  • Python 中的二进制缓冲区

    在Python中你可以使用StringIO https docs python org library struct html用于字符数据的类似文件的缓冲区 内存映射文件 https docs python org library mmap
  • 当玩家触摸屏幕一侧时,如何让 pygame 发出警告?

    我使用 pygame 创建了一个游戏 当玩家触摸屏幕一侧时 我想让 pygame 给出类似 你不能触摸屏幕两侧 的错误 我尝试在互联网上搜索 但没有找到任何好的结果 我想过在屏幕外添加一个方块 当玩家触摸该方块时 它会发出警告 但这花了很长
  • Geopandas 设置几何图形:MultiPolygon“等于 len 键和值”的 ValueError

    我有 2 个带有几何列的地理数据框 我将一些几何图形从 1 个复制到另一个 这对于多边形效果很好 但对于任何 有效 多多边形都会返回 ValueError 请指教如何解决这个问题 我不知道是否 如何 为什么应该更改 MultiPolygon
  • 方程“a + bx = c + dy”的积分解

    在等式中a bx c dy 所有变量都是整数 a b c and d是已知的 我如何找到整体解决方案x and y 如果我的想法是正确的 将会有无限多个解 由最小公倍数分隔b and d 但我只需要一个解决方案 我可以计算其余的 这是一个例
  • ExpectedFailure 被计为错误而不是通过

    我在用着expectedFailure因为有一个我想记录的错误 我现在无法修复 但想将来再回来解决 我的理解expectedFailure是它会将测试计为通过 但在摘要中表示预期失败的数量为 x 类似于它如何处理跳过的 tets 但是 当我
  • Python:尝试检查有效的电话号码

    我正在尝试编写一个接受以下格式的电话号码的程序XXX XXX XXXX并将条目中的任何字母翻译为其相应的数字 现在我有了这个 如果启动不正确 它将允许您重新输入正确的数字 然后它会翻译输入的原始数字 我该如何解决 def main phon
  • 为美国东部以外地区的 Cloudwatch 警报发送短信?

    AWS 似乎没有为美国东部以外的 SNS 主题订阅者提供 SMS 作为协议 我想连接我的 CloudWatch 警报并在发生故障时接收短信 但无法将其发送到 SMS YES 经过一番挖掘后 我能够让它发挥作用 它比仅仅选择一个主题或输入闹钟
  • 设置 torch.gather(...) 调用的结果

    我有一个形状为 n x m 的 2D pytorch 张量 我想使用索引列表来索引第二个维度 可以使用 torch gather 完成 然后然后还设置新值到索引的结果 Example data torch tensor 0 1 2 3 4
  • VSCode:调试配置中的 Python 路径无效

    对 Python 和 VSCode 以及 stackoverflow 非常陌生 直到最近 我已经使用了大约 3 个月 一切都很好 当尝试在调试器中运行任何基本的 Python 程序时 弹出窗口The Python path in your
  • 对输入求 Keras 模型的导数返回全零

    所以我有一个 Keras 模型 我想将模型的梯度应用于其输入 这就是我所做的 import tensorflow as tf from keras models import Sequential from keras layers imp
  • 在 Python 类中动态定义实例字段

    我是 Python 新手 主要从事 Java 编程 我目前正在思考Python中的类是如何实例化的 我明白那个 init 就像Java中的构造函数 然而 有时 python 类没有 init 方法 在这种情况下我假设有一个默认构造函数 就像
  • 您可以在 Python 类型注释中指定方差吗?

    你能发现下面代码中的错误吗 米皮不能 from typing import Dict Any def add items d Dict str Any gt None d foo 5 d Dict str str add items d f
  • Python - 字典和列表相交

    给定以下数据结构 找出这两种数据结构共有的交集键的最有效方法是什么 dict1 2A 3A 4B list1 2A 4B Expected output 2A 4B 如果这也能产生更快的输出 我可以将列表 不是 dict1 组织到任何其他数
  • Python:元类属性有时会覆盖类属性?

    下面代码的结果让我感到困惑 class MyClass type property def a self return 1 class MyObject object metaclass MyClass a 2 print MyObject

随机推荐