# How to write a differential operator in Mathematica

I have a very basic question. I am trying to code the LHS of the differential equation: $$\Big[\frac{1}{\sin\theta}\frac{d}{d\theta}\big(\sin\theta\frac{d}{d\theta}\big)-\frac{m^2}{\sin^2\theta}+n(n+1)\Big]f(\theta)=0$$ into Mathematica. As you can see, in the first term, the differential operator is nested and is to properly act on the function $f(\theta)$.

I tried to use #1, but (lol) I don't know how to use it!

This is what I have:

Sin[θ]^(-1) D[Sin[θ] D[#1, θ], θ] - m^2*#1/Sin[θ]^2 + n (n + 1) #1


What next?!

Also: I don't want to just stick in f[θ] into where I have #1, because eventually I would like to be able to make change of variables!

• Aren't you just missing an & at the end? – Ajasja Mar 2 '13 at 20:48
• One way to go about this is discussed in a notebook available here. See the first example in the section "Some noncommutative algebraic manipulation" – Daniel Lichtblau Mar 2 '13 at 20:54
• @Ajasja Hey! it looks like this is on the right track. I added an & at the end, and all the #1's turned Green. What next? – QuantumDot Mar 2 '13 at 21:08
• You might also be interested in this: Having the derivative be an operator and How to change coordinates of a differential operator? – Jens Mar 2 '13 at 22:25
• @Jens Thank you for these reference! They are very helpful! – QuantumDot Mar 3 '13 at 3:35

Why not just define a Function ? ClearAll[m, n, x, foo]
oper = Function[{f, \[Theta]}, (1/
Sin[\[Theta]] D[Sin[\[Theta]] D[f[\[Theta]], \[Theta]]] -
m^2/Sin[\[Theta]]^2 f[\[Theta]] + n (n + 1) f[\[Theta]])];
oper[foo, x]
foo[x_] := Sin[x]*Cos[x]*x^4
oper[foo, x]

• (+1) That's how I'd do it. – Jens Mar 2 '13 at 22:24
• Except - you don't need the theta in the argument list. In fact, this mixes two alternative approaches: either make the operator always differentiate its argument f with respect to the global variable theta, or do the differentiation with respect to the first slot in f. It's two different things. – Jens Mar 2 '13 at 22:30
• +1, that's the obvious way. but you've mixed up $\theta$ and #1 :) – acl Mar 2 '13 at 22:37
• @Nasser sorry I didn't have time to respond yesterday. I added my own answer to explain what I'd do differently. – Jens Mar 3 '13 at 18:22

Well, I hope you're not pulling my leg here (at 1000 rep one would think you have met pure (anonymous) functions before, or have at least learned where to get help: type & and press F1). But here it goes:

(*Assign the operator in the form of an anonymous function  to oper*)
oper = Sin[\[Theta]]^(-1) D[Sin[\[Theta]] D[#1, \[Theta]], \[Theta]] -
m^2*#1/Sin[\[Theta]]^2 + n (n + 1) #1 &;

(*use it*)
oper [f[\[Theta]]]

• Lol! I am not pulling your leg. My 1000 reputation is coming from asking 1000 idiotic questions. Check my profile. But your response is very helpful! – QuantumDot Mar 3 '13 at 3:34

The differential operator in this question is itself indexed by two variables m and n. This arises after expressing the Laplace operator in spherical coordinates (see the answer by b.gatessucks, but that answer was incomplete because it dropped the $\phi$ dependence too early) and then assuming that the function on which we operate is proportional to $\exp(im\phi)$.

In the other answers that define a function for this operator, the parameters m and n appear as global variables inside the function. That's inconvenient in practice because these indices correspond to physical parameters (the angular momentum components) that are adjustable. So m and n need to appear in the list of variables on which the operator depends.

The next thing on which the operator depends is the independent variable $\theta$ with respect to which we want to take derivatives. It would be reasonable to make this work along the lines of D, which means in particular that it should handle a case like this:

foo[x_] := Sin[θ] Sin[x]
D[foo[θ], θ]


2 Sin[θ] Cos[θ]

This will not work if you define an operator function like this:

myOp = Function[{f, θ}, D[f[θ], θ]];
myOp[foo[θ], θ]


$\sin ^2(\theta )'(\theta )+(2 \sin (\theta )\cos (\theta ))(\theta )$

The problem is that myOp expects to be passed a function instead of an expression. This is not how D works. The purpose of the second argument of D[f,x] is to identify the variable appearing in f that is to be varied. The above method in myOp is basically what Nasser used in his answer. It performs a differentiation with respect to the first slot in the first argument assuming it is a function, and then return the result with the variable $\theta$ inserted in that slot. That's not wrong, it's just an interpretation of differentiation that differs from what D does.

An approach more consistent with D is this definition:

newOp = Function[{f, θ}, D[f, θ]];
newOp[foo[θ], θ]


2 Cos[θ] Sin[θ]

But in this question, the variable θ plays a role that is better grouped together with the indices m and n, and separated from the function f on which the operator acts. Now that I identified how I want to think about the variables on which the differential operator depends, the definition I would use is as follows:

Clear[operNM, n, m, var]
operNM[n_, m_, var_] = Function[{f}, D[Sin[var] D[f, var], var]/Sin[var] -
(m^2 f)/Sin[var]^2 + n (n + 1) f
];


This defines operNM[n, m, var] as an operator to which you then supply an expression that it operates on:

operNM[n,m,θ][foo[θ]]


$\csc (\theta ) \left(4 \sin (\theta ) \cos^2(\theta )-2 \sin ^3(\theta )\right)-m^2+n(n+1) \sin ^2(\theta )$

Having the operator act on expressions instead of expecting its argument to be a function is more practical especially when you deal with functions of several variables. In quantum mechanics, the name of the coordinate in a many-particle wave function is tied to the index of the particle it represents, and function slots are then not a convenient way to keep track of that information. In operNM you can choose the desired differentiation variable independently of where in the expression it appears.

If you are on version 9 you could use the built-in functions (previously available in the VectorAnalysis package). The following might get you close to what you need :

Laplacian[f[theta], {r, theta, phi}, "Spherical"] // Simplify
(* (Cot[theta] Derivative[f][theta] + (f^\[Prime]\[Prime])[theta])/r^2 *)