4

[B]次の[C]ような微分演算子の 2 つの行列が必要です。

B = sympy.Matrix([[ D(x), D(y) ],
                  [ D(y), D(x) ]])

C = sympy.Matrix([[ D(x), D(y) ]])

ans = B * sympy.Matrix([[x*y**2],
                        [x**2*y]])
print ans
[x**2 + y**2]
[      4*x*y]

ans2 = ans * C
print ans2
[2*x, 2*y]
[4*y, 4*x]

これは、次のようなベクトル フィールドのカールの計算にも適用できます。

culr  = sympy.Matrix([[ D(x), D(y), D(z) ]])
field = sympy.Matrix([[ x**2*y, x*y*z, -x**2*y**2 ]])

Sympy を使用してこれを解決するには、次の Python クラスを作成する必要がありました。

import sympy

class D( sympy.Derivative ):
    def __init__( self, var ):
        super( D, self ).__init__()
        self.var = var

    def __mul__(self, other):
        return sympy.diff( other, self.var )

微分演算子の行列が左側で乗算されている場合、このクラスだけで解決します。ここでdiffは、微分する関数がわかっている場合にのみ実行します。

微分演算子の行列が右側で乗算されている場合の回避策__mul__として、コア クラスのメソッドをExpr次のように変更する必要がありました。

class Expr(Basic, EvalfMixin):
    # ...
    def __mul__(self, other):
        import sympy
        if other.__class__.__name__ == 'D':
            return sympy.diff( self, other.var )
        else:
            return Mul(self, other)
    #...

これはかなりうまく機能しますが、これを処理するためのより良いネイティブ ソリューションが Sympy にあるはずです。誰がそれが何であるか知っていますか?

4

4 に答える 4

5

このソリューションは、他の回答とここからのヒントを適用します。演算子は次のDように定義できます。

  • 左から乗算された場合のみ考慮されるため、何もD(t)*2*t**3 = 6*t**22*t**3*D(t)ません
  • で使用されるすべての式と記号はDis_commutative = False
  • を使用して、指定された式のコンテキストで評価されますevaluateExpr()
    • これは式に沿って右から左に移動し、演算子を見つけて* を対応する右側の部分にD適用しますmydiff()

*:のように、より高い順序を作成できるようにするmydiff代わりに使用されます。diffDmydiff(D(t), t) = D(t,t)

現在のソリューションでは が実際に微分作業を行っているため、diff内部__mul__()は参照用にのみ保持されています。python mudule が作成され、 として保存されました。DevaluateExpr()d.py

import sympy
from sympy.core.decorators import call_highest_priority
from sympy import Expr, Matrix, Mul, Add, diff
from sympy.core.numbers import Zero

class D(Expr):
    _op_priority = 11.
    is_commutative = False
    def __init__(self, *variables, **assumptions):
        super(D, self).__init__()
        self.evaluate = False
        self.variables = variables

    def __repr__(self):
        return 'D%s' % str(self.variables)

    def __str__(self):
        return self.__repr__()

    @call_highest_priority('__mul__')
    def __rmul__(self, other):
        return Mul(other, self)

    @call_highest_priority('__rmul__')
    def __mul__(self, other):
        if isinstance(other, D):
            variables = self.variables + other.variables
            return D(*variables)
        if isinstance(other, Matrix):
            other_copy = other.copy()
            for i, elem in enumerate(other):
                other_copy[i] = self * elem
            return other_copy

        if self.evaluate:
            return diff(other, *self.variables)
        else:
            return Mul(self, other)

    def __pow__(self, other):
        variables = self.variables
        for i in range(other-1):
            variables += self.variables
        return D(*variables)

def mydiff(expr, *variables):
    if isinstance(expr, D):
        expr.variables += variables
        return D(*expr.variables)
    if isinstance(expr, Matrix):
        expr_copy = expr.copy()
        for i, elem in enumerate(expr):
            expr_copy[i] = diff(elem, *variables)
        return expr_copy
    return diff(expr, *variables)

def evaluateMul(expr):
    end = 0
    if expr.args:
        if isinstance(expr.args[-1], D):
            if len(expr.args[:-1])==1:
                cte = expr.args[0]
                return Zero()
            end = -1
    for i in range(len(expr.args)-1+end, -1, -1):
        arg = expr.args[i]
        if isinstance(arg, Add):
            arg = evaluateAdd(arg)
        if isinstance(arg, Mul):
            arg = evaluateMul(arg)
        if isinstance(arg, D):
            left = Mul(*expr.args[:i])
            right = Mul(*expr.args[i+1:])
            right = mydiff(right, *arg.variables)
            ans = left * right
            return evaluateMul(ans)
    return expr

def evaluateAdd(expr):
    newargs = []
    for arg in expr.args:
        if isinstance(arg, Mul):
            arg = evaluateMul(arg)
        if isinstance(arg, Add):
            arg = evaluateAdd(arg)
        if isinstance(arg, D):
            arg = Zero()
        newargs.append(arg)
    return Add(*newargs)

#courtesy: https://stackoverflow.com/a/48291478/1429450
def disableNonCommutivity(expr):
    replacements = {s: sympy.Dummy(s.name) for s in expr.free_symbols}
    return expr.xreplace(replacements)

def evaluateExpr(expr):
    if isinstance(expr, Matrix):
        for i, elem in enumerate(expr):
            elem = elem.expand()
            expr[i] = evaluateExpr(elem)
        return disableNonCommutivity(expr)
    expr = expr.expand()
    if isinstance(expr, Mul):
        expr = evaluateMul(expr)
    elif isinstance(expr, Add):
        expr = evaluateAdd(expr)
    elif isinstance(expr, D):
        expr = Zero()
    return disableNonCommutivity(expr)

例 1: ベクトル フィールドのカール。commutative=False変数の順序がMul().args結果に影響するため、変数を定義することが重要であることに注意してください。この他の質問を参照してください。

from d import D, evaluateExpr
from sympy import Matrix
sympy.var('x', commutative=False)
sympy.var('y', commutative=False)
sympy.var('z', commutative=False)
curl  = Matrix( [[ D(x), D(y), D(z) ]] )
field = Matrix( [[ x**2*y, x*y*z, -x**2*y**2 ]] )       
evaluateExpr( curl.cross( field ) )
# [-x*y - 2*x**2*y, 2*x*y**2, -x**2 + y*z]

例 2: 構造解析で使用される典型的なリッツ近似。

from d import D, evaluateExpr
from sympy import sin, cos, Matrix
sin.is_commutative = False
cos.is_commutative = False
g1 = []
g2 = []
g3 = []
sympy.var('x', commutative=False)
sympy.var('t', commutative=False)
sympy.var('r', commutative=False)
sympy.var('A', commutative=False)
m=5
n=5
for j in xrange(1,n+1):
    for i in xrange(1,m+1):
        g1 += [sin(i*x)*sin(j*t),                 0,                 0]
        g2 += [                0, cos(i*x)*sin(j*t),                 0]
        g3 += [                0,                 0, sin(i*x)*cos(j*t)]
g = Matrix( [g1, g2, g3] )

B = Matrix(\
    [[     D(x),        0,        0],
     [    1/r*A,        0,        0],
     [ 1/r*D(t),        0,        0],
     [        0,     D(x),        0],
     [        0,    1/r*A, 1/r*D(t)],
     [        0, 1/r*D(t), D(x)-1/x],
     [        0,        0,        1],
     [        0,        1,        0]])

ans = evaluateExpr(B*g)

大きな式をすばやくチェックするためのprint_to_file()関数が作成されました。

import sympy
import subprocess
def print_to_file( guy, append=False ):
    flag = 'w'
    if append: flag = 'a'
    outfile = open(r'print.txt', flag)
    outfile.write('\n')
    outfile.write( sympy.pretty(guy, wrap_line=False) )
    outfile.write('\n')
    outfile.close()
    subprocess.Popen( [r'notepad.exe', r'print.txt'] )

print_to_file( B*g )
print_to_file( ans, append=True )
于 2013-04-26T12:57:54.797 に答える
3

微分演算子は SymPy のコアには存在せず、存在したとしても「演算子の適用」ではなく「演算子による乗算」は、SymPy でサポートされていない表記の乱用です。

[1] もう 1 つの問題は、SymPy 式は のサブクラスからしか構築できないため、 として入力すると単純にエラーが発生するsympy.Basic可能性が高いことです。これが失敗する理由です。構築できません。class Dsympy_expr+D(z)(expression*D(z)) * (another_expr)(expression*D(z))

さらに、 の引数がD単一Symbolでない場合、この演算子に何を期待するかが明確ではありません。

最後に、diff(f(x), x)(wherefはシンボリックな未知の関数です) は、あなたが観察したように、評価されていない式を返しますf。後で、expr.subs(f(x), sin(x))導関数を代入すると評価されます (最悪の場合、 を呼び出す必要があるかもしれませんexpr.doit())。

[2] あなたの問題に対するエレガント短い解決策はありません。あなたの問題を解決するために私が提案する方法は、の__mul__メソッドをオーバーライドすることですExpr:式ツリーを乗算するだけでなく、左の式ツリーにインスタンスが含まれているかどうかを確認し、Dそれらを適用します。新しいオブジェクトを追加したい場合、明らかにこれはスケーリングしません。これは、sympy の設計に関する長年の既知の問題です。

編集: [1] は、単に を含む式の作成を許可するために必要ですD。[2] は、1 つ以上の何かを含む式が機能するために必要ですD

于 2013-03-17T19:39:43.147 に答える
1

正しい乗算を機能させたい場合は、からサブクラス化する必要がありますobject。それはx*DにフォールバックしD.__rmul__ます。ただし、演​​算子が右から適用されることはないため、これが優先度が高いとは想像できません。

于 2013-03-18T17:03:35.263 に答える
1

常に自動的に機能するオペレーターを作成することは、現時点では不可能です。本当に完全に機能するには、 http://code.google.com/p/sympy/issues/detail?id=1941が必要です。https://github.com/sympy/sympy/wiki/Canonicalizationも参照してください(そのページは自由に編集してください)。

ただし、そのstackoverflowの質問からのアイデアを使用して、ほとんどの場合機能するクラスを作成できます。それが処理されない場合は、式を通過し、演算子が適用されていない場所に演算子を適用する単純な関数を記述しますまだ適用されていません。

ところで、「乗算」としての微分演算子で考慮すべきことの 1 つは、それが非結合であることです。つまり、(D*f)*g=g*Dfに対してD*(f*g)=g*Df + f*Dgです。したがって、式全体ではなく一部を「食べる」ことがないように注意する必要があります。たとえば、このためにD*2*x与えるでしょう0。SymPy はどこでも乗算が連想的であると想定しているため、ある時点でそれが正しく行われない可能性があります。

それが問題になる場合は、自動アプリケーションをダンプし、それを通過して適用する関数で作業することをお勧めします (上記で述べたように、とにかく必要になります)。

于 2013-03-18T17:15:26.377 に答える