Symbolic rewriting, expanding and collecting

Specific Domains
, ,
1

Hi!

I need to revise a lot of mundane analytical calculations which I guess could be automated. In its essence, I have three sets of symbols: operators E, X, C, S \in \rm{Op}, sigma matrices \eta_x, \eta_y, \eta_z \in SU(2) and a set of constants a,c,s \in \rm{Re}. These symbols have the following rules:

  • \eta_x \eta_y = i \eta_z, \eta_y \eta_z = i \eta_x, \eta_z \eta_x = i\eta_y, \eta_y \eta_x =- i \eta_z, \eta_z \eta_y = -i \eta_x, \eta_x \eta_z = -i\eta_y
  • The symbols of Op do not commute with themselves
  • Symbols from Re commute with themselves
  • The symbols from different sets commute with each other

With these rules, I need to expand multiplications such as:

(E \eta_z - a X \eta_x) (C - \frac{i\eta_y}{2}\ S)(c \eta_x - s \eta_y ) (C + \frac{i \eta_y}{2} S) = i c \eta_y E (C + \frac{i \eta_y}{2} S) + i s \eta_x E - a c X (C + \frac{i \eta_y}{2} S) + i s \eta_z X

which is a simple but rather lengthy calculation which is hard to make manageable. Are there any Julia package or any other CAS which could help me with such error-prone tasks?

2

SymbolicUtils.jl would be the way to go with this. Define a ruleset and it’ll apply it.

3 Likes
3

SymbolicUtils looks cool. However, I cannot understand how to use it for expressions that are not with predefined types Real, Complex and Int. Is it possible to define nonnumber symbols or types?

4

It’s literally just a dispatching rule system, so you just create a different set of rules.

Nothing specifies Number.

1 Like
5

It is quite hard to get into the meaning of these rules. I guess my biggest obstacle is what is the meaning of ~? In the Github README.md, I see that the rules are defined with single ~. In the code, I see ~~x to be used most often.

I see it is also possible to define rules without ~ for rules applying to particular symbols. I guess it would be the right place do define rules for \eta_{x,y,z}. However, I am puzzled why for instance if I define a rule

 r = @rule sinh(im * y) => sin(y)

and apply it to the example

simplify(cos(y)^2 + sinh(im*y)^2, RuleSet([r]))

I get (cos(y))^2 + (sin(y))^2 instead of 1 if the rule would be defined with ~x.

Also, what is the difference between @acrule and @rule?

6

Commutativity IIRC. It’s all explained in a bit of detail:

https://juliasymbolics.github.io/SymbolicUtils.jl/api/#rewriters

2 Likes
7

Open an issue for any documentation you feel is missing so @shashi will see it.

1 Like
8

Looking into the library and trying to implement Pauli matrix algebra, I realised that rule before applied seems to sort its arguments under multiplications. For instance, I can not get a rule to apply on the following code:

abstract type SU2 end
@syms ηx::SU2 ηy::SU2 ηz::SU2
(@rule ηy*ηx => -im*ηz)(ηy*ηx)==nothing

which returns nothing. Instead, it applies the wrong rule:

(@rule ηx*ηy => im*ηz)(ηy*ηx)

Ahh, in fact, it sorts all multiplications. Is there any way to create noncommuting symbols under multiplication?

9

open an issue.

10

Replied with the solution here Symbols gets sorted even for user types · Issue #206 · JuliaSymbolics/SymbolicUtils.jl · GitHub

1 Like
11

So I tried to work around the issue by defining a \times method for noncommutative terms. For a moment, I try to implement only Pauli operator algebra and forget at the moment about the set of operators \rm{Op}. I can get substitution rules to work in some situation, whereas in others, it does fail. I came up with the following code to define simplification rules for operators \eta_{x,y,z}:

using SymbolicUtils
using SymbolicUtils: Sym, Mul, Term, Symbolic, Add

abstract type SU2 end

@syms ηx::SU2 ηy::SU2 ηz::SU2 α::Real β::Real γ::Real

×(x::Symbolic, y::Symbolic) = Term(*, [x, y])

function ×(x::Symbolic, y::Term) 
	if y.f==*
		Term(*, [x, y.arguments...])
	else
		Term(*, [x, y])
	end
end

function ×(x::Term, y::Symbolic) 
	if x.f==*
		Term(*, [x.arguments..., y])
	else
		Term(*, [x, y])
	end
end

×(x::Add, y::Symbolic) = Term(*, [x, y])

srules = [
    @rule ηx*ηy => im*ηz
    @rule ηy*ηz => im*ηx
    @rule ηz*ηx => im*ηy
    @rule ηy*ηx => -im*ηz
    @rule ηz*ηy => -im*ηx
    @rule ηx*ηz => -im*ηy
	
	@rule ηx^2 => 1
    @rule ηy^2 => 1
    @rule ηz^2 => 1
	
    @rule ηx*ηy*~x => im*ηz*~x
    @rule ηy*ηz*~x => im*ηx*~x
    @rule ηz*ηx*~x => im*ηy*~x
    @rule ηy*ηx*~x => -im*ηz*~x
    @rule ηz*ηy*~x => -im*ηx*~x
    @rule ηx*ηz*~x => -im*ηy*~x
    
    @rule ~x*ηx*ηy => ~x*im*ηz
    @rule ~x*ηy*ηz => ~x*im*ηx
    @rule ~x*ηz*ηx => ~x*im*ηy
    @rule ~x*ηy*ηx => ~x * (-im)*ηz
    @rule ~x*ηz*ηy => ~x * -im*ηx
    @rule ~x*ηx*ηz => ~x * -im*ηy
	
	@rule (~x + ~y)*~z => ~x*~z + ~y*~z
	@rule ~z*(~x + ~y) => ~z*~x + ~z*~y
    
]

Note I repeated the rules with ~x as otherwise these rules were not substituted.

Here are some tests I used to evaluate the rules:

simplify(ηx×ηy×ηz + ηz×ηx, RuleSet(srules)) # Success

simplify(ηx×ηy×ηz×ηy, RuleSet(srules)) # Failure

simplify((ηz×ηy + ηz)×ηz, RuleSet(srules)) # Success/Failure

simplify(α×ηx×ηy×ηz + β×ηz×ηx, RuleSet(srules)) # Works only for the second term

simplify(ηx×ηy×ηz×α + ηz×ηx×β, RuleSet(srules)) # Works only for the second term

simplify((α×ηx + β×ηz×ηx)×(γ×ηy + ηz), RuleSet(srules)) # Failure
12
  • There is a missing x(::Term, ::Term) method which may create ambiguity.
  • *(x::Symbolic{SU2}, y::Symbolic{SU2}) = Term(*, [x, y]) might be necessary because you are using * in the RHS of your rules. You may consider this or use × even in the rules.