I have the functional:

<https://lh5.googleusercontent.com/-YyVcmvfWlB8/UnIq_SuBdjI/AAAAAAAABYE/NkiWxG35LS8/s1600/fig1.png>

Where A is a function of "v". The non-linear system of equations necessary 
to find "v" is obtained doing:

<https://lh5.googleusercontent.com/-SElccyNi1xo/UnIrNKzuKxI/AAAAAAAABYM/xd6tHk0SJT0/s1600/fig2.png>

According to the differentiation rules, this system could be written in 
matrix form as:

<https://lh6.googleusercontent.com/-Qm0cHaMnEAA/UnIrYLIfs5I/AAAAAAAABYU/ocDPAfGA4ic/s1600/fig3.png>

So I've implemented in Sympy, with the help of NumPy "ndarray" the 
differentiation about a vector. In this case the shape of (\partial A / 
\partial v is (n X n X n), where "n" is the length of vector "v".

The function is simple:
def vdiff(x, vector):
    x = np.array(x)
    shape = x.shape
    ans = [np.array([e.diff(vi) for e in x.ravel()]) for vi in vector]
    ans = [a.reshape(shape) for a in ans]
    return np.array(ans).swapaxes(0, 1)

 

And I've tested with two different cases:

def test_00():
    from sympy.abc import a, b, c, x
    A11 = x*a*c**2
    A12 = x**2*a*b*c
    A13 = x**2*a**3*b**5
    A21 = x**3*a**2*b*c
    A22 = x**4*a*b**2*c**5
    A23 = 5*x**4*a*b**2*c
    A31 = x**4*a*b**2*c**4
    A32 = 4*x**4*a*b**2*c**2
    A33 = 4*x**4*a**5*b**2*c
    A = np.array([[A11, A12, A13],
                  [A21, A22, A23],
                  [A31, A32, A33]])
    v = np.array([a, b, c])
    F = (v.T.dot(A)).dot(v)
    Av = vdiff(A, v)
    p1 = v.dot(A)
    p2 = A.dot(v)
    p3 = v.dot(Av.dot(v))
    new = p1 + p2 + p3
    ref = np.array([F.diff(a), F.diff(b), F.diff(c)])
    print sympy.simplify(Matrix(ref-new))

def test_01():
    from sympy.abc import a, b, c, x
    sympy.var('c1, c2')
    A11 = x*a*c**2
    A12 = x**2*a*b*c
    A13 = x**2*a**3*b**5
    A21 = x**3*a**2*b*c
    A22 = x**4*a*b**2*c**5
    A23 = 5*x**4*a*b**2*c
    A = np.array([[A11, A12, A13],
                  [A21, A22, A23]])
    v = np.array([a, b, c])
    cc = np.array([c1, c2])
    F = cc.dot(A.dot(v))
    ref = np.array([F.diff(a), F.diff(b), F.diff(c)])
    Av = vdiff(A, v)
    p1 = cc.dot(A)
    p2 = cc.dot(Av.dot(v))
    new = p1 + p2
    print sympy.simplify(Matrix(ref-new))


I'd like to share here and get some feedback in case someone has a more 
straightforward way to implement this...


Thank you!

Saullo




-- 
You received this message because you are subscribed to the Google Groups 
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sympy+unsubscr...@googlegroups.com.
To post to this group, send email to sympy@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.

Reply via email to