Recursive function for Integration by parts algorithm

Question

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 get the answer 0, instead of precisely the missing terms.

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)

Answer 1

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]]