# Recursive function for Integration by parts algorithm

I need to write a function that performs repeatedly integration by parts on an integral, in views of transferring the derivatives from one function to the other. For example, in the case where the second function is a second order derivative, I want to obtain automatically the right hand side below, where derivatives have been transferred to $u$:

$$\iint u(r,\theta)\,\partial_{rr}\psi(r,\theta)\,d^2 r=u\psi-2\int\psi\,\partial_r u\,dr+\iint\psi\,\partial_{rr}u\,d^2 r$$

To this end, I have written the following Mathematica code, where IPs is a function defined recursively, given below:

IPs[1, u_, v_, r] := u*Integrate[v, r] - Integrate[D[u, r]*Integrate[v, r], r];
IPs[n_Integer, u_, v_, r] /; n > 1 :=
IPs[n - 1, u, Integrate[v, r], r] - Integrate[IPs[n - 1, D[u, r], Integrate[v, r], r], r]


I test this function by defining two functions which depend formally on two variable, but really depend on the first variable only.

f[x_, y_] := u[x]
g[x_, y_] := D[D[ψ[x], x], x]
IPs[2, f[r, θ], g[r, θ], r]


At this point, I should get the $\TeX$ formula above. Instead, I obtain:

$$\iint u(r,\theta)\,\partial_{rr}\psi(r,\theta)\,d^2 r=u\psi-\int\psi\,\partial_r u\,dr$$

Half of the terms are missing. More simply, if I invoke

Integrate[IPs[1, f[r, θ], g[r, θ], r], r]


I am relatively green in Mathematica. My best guess is that, instead of returning a function, IPs returns a different kind of object that doesn't interact well with Mathematica's Integrate.

What am I doing wrong here?

NB: I have considered using nested functions, but that does not seem helpful for the task I want to accomplish.

NB2: Using SymPy, the following Python code works (at the cost of a very poor legibility of the output, and at a very slow speed).

from sympy import *

x=symbols('x')
y=symbols('y')

def f(x,y):
u = Function('u')
return u(x)

def g(x,y):
p = Function('p')
return diff(diff(diff(p(x),x),x),x)

def IPs(n,a,b,r):
if n==1:
ans = a * integrate(b,r) - integrate(diff(a,r)*integrate(b,r),r)
return ans
else:
ans = IPs(n-1,a,integrate(b,r),r) - integrate(IPs(n-1,diff(a,r),integrate(b,r),r),r)
return ans

dummy=IPs(3,f(x,y),g(x,y),x)
print(dummy)

• The problem can be boiled down to the following: Integrate[-Integrate[\[Psi][r] u''[r], r] + \[Psi][r]*u'[r], r] gives 0. Looks like a bug. – xzczd Oct 6 '17 at 4:38
• This seems to be a known bug: mathematica.stackexchange.com/q/37109/1871 – xzczd Oct 6 '17 at 4:45

As mentioned above, you seem to encounter a bug of Integrate. In your case,

Integrate[-Integrate[ψ[r] u''[r], r] + ψ[r]*u'[r], r]


mistakenly evaluates to 0. However, even if this bug doesn't exist, your function won't work well, because Integrate doesn't automatically apply linearity to integrand in this case. (For example, Integrate[f'[x] + f[x], x] won't evaluate to f[x] + Integrate[f[x], x]. )

To circumvent the problem, we need to implement the algorithm in a way doesn't rely on Integrate that much. The following is my implementation:

Clear@int
int[deri_ rest_, var_, varrest___, deri_] /; ! FreeQ[deri, Derivative] :=
With[{intderi = Integrate[deri, var]},
int[intderi rest, varrest, intderi] - int[D[rest, var] intderi, var, varrest, intderi]]
int[expr_, deri_] := expr
Format@int[expr_, var__, deri_] := Integrate[expr, var]

int[u[r, t] D[ψ[r, t], r, r], r, r, D[ψ[r, t], r, r]]


• Interesting bug. I am surprised that mathematica can fail so spectacularly on a high school problem. But fine, bugs exist. Thanks for your suggestion. To be honest I have a hard time following the flow here, but it does work to eliminate second order derivatives. Could your function be invoked to get rid of third order derivatives, as in my sympy example above? (Or any arbitrary order, which is really the ultimate goal) Trying int[u[r, t] D[[Psi][r, t], r, r, r], r, r, D[[Psi][r, t], r, r, r]] still returns first order derivatives of 'Psi' – BenBoulderite Oct 6 '17 at 15:54