What is sympy (python) equivalent differentiation method in Julia?

General Usage
1

i have been struggling to find out what is the equivalent sympy method in Julia. actually i made finite element method solver (numerical method to solve partial differential equation, parabolic, elliptical and hyperbolic in 2D) and it is working fine in python. but, i am facing the speed issue, as this solver contain very large number of string operation i mean string expression. i will mention brief steps involved in solver as follows:

  1. make string shape functions for each rectangular element in array
  2. making functional of above string shape function as per Residual Galerkin method
  3. find jacobian of given physical element in standard rectangular element
  4. taking differentiation (symbolically and not numerically)–here i am facing problem in julia
  5. taking integration of each functional element wise (numerically)
  6. obtaining Stiff, Mass, and Load vectors
  7. solving matrices equation using linear aljebra
  8. post processing
    end
    i want to make same solver in Julia to increase the speed. please help to make solver in Julia.

i am sending some snippet from solver.

def vectorshapes(p):
“”" Return vector shape functions of edges for order p, first all horizontal edges then vertical edges
-return list of touple [(ux,uy)…]
“”"
s,t=symbols(‘s t’)
shape=masterRectShapes(p)# gives Lagrange shape functions in string format containg s and t variables.
HorizontalEdges =[(i+(p+1)j,i+(p+1)(j+1)) for i in range(p+1) for j in range(p)]
VerticalEdges =[((p+1)i+j,(p+1)i+(j+1)) for i in range(p+1) for j in range(p)]
HorizontalEdges.extend(VerticalEdges)
del VerticalEdges
edges=HorizontalEdges
edges=[(f’(2/{pp})(({shape[i]}({diff(shape[j],s)}))-({shape[j]}({diff(shape[i],s)})))‘,f’(2/{pp})(({shape[i]}({diff(shape[j],t)}))-({shape[j]}({diff(shape[i],t)})))') for i,j in edges]
return edges
here i have to differentiate shape fuction wrt s for x components and wrt t for y components

2

Just two general remarks:

  • if you post code snippets, please enclose them in triple backticks (```) for better readability
  • there is not only numeric and symbolic differentiation, but also automated differentiation which is more common in Julia code. Is that something you could use?

You could also look at: Derivatives and Differentials · Symbolics.jl , but this library is not as mature as SymPy …

3 Likes
3

thanks for your response, the finite element method is used for arbitrary domain and not for regular domain. i can not use simple integration and derivative technique as used for regular domain. my problem is intermediate step only, that is, i have to use this symbolic differentiation technique for this function only. my objective is to find for example, dx^2=2*x in in symbolic form, so in next step i want to use these result for numerical integration using Gauss Quadrature rules.

4

I want derivatives like f(“x^2y", x)=2xy and f("x^2y”, y)=x^2 output is string or symbol
(partial derivative wrt x and y), any function in julia like this?

5

SymPy works with symbolic expressions, not strings. You maybe can use sympify, as in diff(sympify("x^3*y"), x), though that will be an issue if you have assumptions on x:

julia> @syms x y
(x, y)

julia> diff(sympify("x^3*y"), x)
   2  
3⋅x ⋅y

julia> @syms x::real
(x,)

julia> diff(sympify("x^3*y"), x)
0
6

Don’t do this. Represent functions using functions (including anonymous functions), not strings.

Functions are fast, composable, flexible (exposing the full Julia language), and differentiable using any of the many Julia AD packages. Strings (as representations of mathematical expressions) are the opposite of these things.

(Check out Gridap.jl, which is a full-featured Julia finite-element package.)

2 Likes
7

thanks, i welcome your suggestion but, if i had not made strings shape function then my solver would have particular for single differential equation and particular boundary and initial condition(in time dependent). i made this solver in such way that, just you have to modify functional in in string form, apply Dirichlet or Neuman or robin condition or all in (string form input) and input vertices of given domain. solution is ready for analysis!. btw, without mesh!
you can think like JSON object work so that others don’t need to do any programming part, he has to just modify functional and boundary condition in string form.

8

thanks, i will try it. will it have speed of Julia? or is it just pycall?

9

You can have function as inputs of another function. Much more efficient in Julia than carrying a string around. So if the argument of your function is "3*x^2", you could just pass x-> 3*x^2
With respect to differentiation, with automatic differentiation you get an exact evaluation of the derivative, its not numerical. The only thing that you need is the function in Julia.

10

and what could happens if user changed the functional like "3x^2(yt^2)(other gradient, divergence, curl etc)?, i am talking about functional and not function.

11

for example,

solver.setLoad('0')
# solver.setRobinCo({1:'1.5'}) # it is dictionary {elementNumber,coefficient value}
# solver.setStiffCo({3:'0.4'})
# solver.setMassCo({2:'2.55'})
startTime=time.time()
totalElements,npe,cords,lg,dirich,neum, robin=mesh.data('eigen.txt')
solver.setFunctional(f"delfx('$i',s,t,j)*delfx('$j',s,t,j)+delfy('$i',s,t,j)*delfy('$j',s,t,j)",p)#$i and $j are shape functions
original_order=int(math.sqrt(npe)-1) # calculation of original order from nodes per element.
elementsCords, npe=solver.getHigherNodes(cords,p,original_order)#[[(),(),()],[(),(),()]]
stiffMatrices,loads=solver.getStiffsLoads(elementsCords,p)# each element cordinates, order 
massMatrices=solver.getMassMatrices(elementsCords,p)#if  you want to mass matrix open this
# loads=solver.setNeumannBC(elementsCords,loads,neum,p)
# stiffMatrices,loads=solver.setRobinBC(elementsCords,stiffMatrices,loads,robin,p,timeInterval)#if problmes occurs there may be in robinBC
globalNumbers,globalCords,totalNodes=solver.getGlobalData(elementsCords)
globalStiff,globalMass,globalLoad=solver.assemble(stiffMatrices,loads,globalNumbers,totalNodes,p,massMatrices)
GK,GM,GL=solver.setDirichletBC(elementsCords,globalStiff,globalLoad,globalNumbers,dirich,p,None,globalMass)
#imp: eigen vectors are each  columns of given 2d array. 
eigv,eigvec = eigh(GK,GM)```

in above program you can make your own functional from any differential equation, and we can change boundary condition in text file (in string form)
12

Depending on that type of functions that you are using on the boundary is what going to be better.

If f is a really simple function, go for symbolic differentiation. But it f is a more complicated function (multiple compositions for instance), just use automatic differentiation.Is much more exact and faster.

With respect to the functional, the functional itself can be an argument of your solver. Just define with nice symbols. Remember that in Julia you can use mathematical symbols.

Also, please read about automatic differentiation using dual numbers, super simple idea but works great in this context.

13

If speed is a problem, then sympy won’t help you, since it’s calling Python. Symbolics.jl/SciML-ecosystem should be very fast for this.

But what are the limitations of automatic differentiation for your use case? Your code example was a bit massive and confusing.

14

I will give simple steps, where my string functional start and how they end at numerical value

  1. suppose i want to solve Laplace equation in 2D, -∇.(c∇u)=0,in Ω
  2. i have decomposed into functional format in string form like
    (f"delfx(‘$i’,s,t,j)*delfx(‘$j’,s,t,j)+delfy(‘$i’,s,t,j)*delfy(‘$j’,s,t,j)",p)
  3. here $i and $j are Lagrange basis function (polynomial contain s and t natural coordinates, in string format)
  4. here coefficient c from original pde, may have function of (x,y,time) then i have to convert these coefficient in to natural coordinates (s,t) then i have to multiply these coefficient to above functional(again they are in string form)
  5. now i have to take numerical integration using Gauss Quadratue rule, for that i have to pass whole functional (string form )to integrate function, here i have to parse and evaluate each term from string at each points in domain with calling jacobian for each points, this jacobian required for integration as well as differentiation.
  6. here i will get tensor of each element and they i can be solved using any linear Aljebra pack.

my approach is to make single program for all differential equation and i did it in python(no problem up to 5 th order, for fifth order program take 15s to… depend on number of element,). in this approach i can not make independent function for particular job for example, if i have to integrate above functional i can not differentiate first and then integrate, i have to call one by one derivative function in integration function with jacobian for each points, this is happing because due to generalization for all pde.
as you mentioned i will read completely about AD Symbolics.jl lmethods in Julia then i will shift from python to Julia.

15

is there any package or codes that solve any differential equation like elliptical, parabolic (time dependent and hyperbolic with three boundary condition and for higher order ) for arbitrary cross section without mesh like gmsh etc in 20-30 lines just by modifying functional, boundary condition vertices of domain and boundary condition? i have made this program for higher accuracy of solution for by increasing the order of element know as hp_finite element method

16

No, I meant that it wasn’t clear what the central, critical part was, where AD could replace string manipulation for derivation. It doesn’t need all the parts to explore just that.

Sounds good. Just note that AD and symbolic differentiation with Symbolics.jl are two different approaches.

17

You can do this without passing functions as strings, as we’ve been trying to tell you. For example, see the Gridap tutorial.

18

For your consideration