Poor performance when evaluating double integral of hat functions

Question

I am trying to evaluate a double integral but am having performance issues.

First a bit of math

I am using a Galerkin like scheme to discretize operators in an equation I am working with. The problem I am trying to solve is this:

I start with hat functions in x and y,

$ \Psi_j(z) \left\{ \begin{array}{ll} \frac{1}{\Delta z}(z-z_{j-1}) & , z < z_j\\ -\frac{1}{\Delta z}(z-z_{j-1}) & , z >z_j \\ 0 & , \left\vert z - z_j \right\vert > \Delta z \\ \end{array} \right. , $

$ \Psi_i(y) \left\{ \begin{array}{ll} \frac{1}{\Delta y}(y-y_{i-1}) & , y < y_i\\ -\frac{1}{\Delta y}(y-y_{i-1}) & , y >y_i \\ 0 & , \left\vert y - y_i \right\vert > \Delta y \\ \end{array} \right. . $

where $i$ and $j$ correspond to a node on our grid. I then define an arbitrary function and the variable of interest (pressure in my case) as,

$ \Theta(y,z) \cong \sum\limits_i \sum\limits_j \Theta_{ij} \Psi_i(y)\Psi_j(z) , $

$ p(y,z) \cong \sum\limits_i \sum\limits_j p_{ij} \Psi_i(y)\Psi_j(z) . $

The arbitrary function corresponds to the coefficient of the operator we are approximating. The integral I am trying to solve is,

$ \Theta p |_{z = z_i} \cong \frac{\int \int \Psi_i \Psi_j \Theta p dz dy}{\int \int \Psi_i \Psi_j dz dy} $

Now the Code

I am building off of what I had done earlier.

Here are the piecewise functions:

ψzt[z_, c_] := Piecewise[{{(z - c)/Δz + 1, z <= c}, 
                          {-(z - c)/Δz + 1, z > c}}];
ψz[z_, c_] := Piecewise[{ {ψzt[z, c], ψzt[z, c] > 0}, 
                          {0, ψzt[z, c] <= 0}}];
ψyt[y_, c_] := Piecewise[{{(y - c)/Δy + 1, y <= c}, 
                          {-(y - c)/Δy + 1, y > c}}];
ψy[y_, c_] := Piecewise[{ {ψyt[y, c], ψyt[y, c] > 0}, 
                          {0, ψyt[y, c] <= 0}}];

Here are the functions I have defined (I just expand the double sum):

θ [y_, z_] := 
Subscript[θ, i - 1, j - 1] ψy[y, Subscript[Y, i - 1]] ψz[z,Subscript[Z, i - 1]] +               
Subscript[θ, i - 1,j] ψy[y, Subscript[Y, i]] ψz[z, Subscript[Z, i - 1]] +  
Subscript[θ, i - 1,j + 1] ψy[y, Subscript[Y, i + 1]] ψz[z, Subscript[Z, i - 1]] +  
Subscript[θ, i,j - 1] ψy[y, Subscript[Y, i - 1]] ψz[z, Subscript[Z,i]] +
Subscript[θ, i,j] ψy[y, Subscript[Y, i]] ψz[z, Subscript[Z, i]] + 
Subscript[θ, i, j + 1] ψy[y, Subscript[Y, i + 1]] ψz[z, Subscript[Z,i]] +  
Subscript[θ, i + 1,j - 1] ψy[y, Subscript[Y, i + 1]] ψz[z, Subscript[Z, i - 1]] +  
Subscript[θ, i + 1,j] ψy[y, Subscript[Y, i]] ψz[z, Subscript[Z, i + 1]] + 
Subscript[θ, i + 1,j + 1] ψy[y, Subscript[Y, i + 1]] ψz[z,Subscript[Z,i + 1]]; 


μ[y_, z_] := 
Subscript[μ, i - 1, j - 1] ψy[y, Subscript[Y, i - 1]] ψz[z, Subscript[Z, i - 1]] +  
Subscript[μ, i - 1, j] ψy[y, Subscript[Y, i]] ψz[z, Subscript[Z, i - 1]] + 
Subscript[μ, i - 1,j + 1] ψy[y, Subscript[Y, i + 1]] ψz[z, Subscript[Z,i - 1]] +  
Subscript[μ, i,j - 1] ψy[y, Subscript[Y, i - 1]] ψz[z, Subscript[Z, i]] +
Subscript[μ, i, j] ψy[y, Subscript[Y, i]] ψz[z, Subscript[Z, i]] +
Subscript[μ, i, j + 1] ψy[y, Subscript[Y, i + 1]] ψz[z, Subscript[Z, i]] +
Subscript[μ, i + 1, j - 1] ψy[y, Subscript[Y, i + 1]] ψz[z, Subscript[Z,i - 1]] + 
Subscript[μ, i + 1,j] ψy[y, Subscript[Y, i]] ψz[z, Subscript[Z, i + 1]] + 
Subscript[μ, i + 1,j + 1] ψy[y, Subscript[Y, i + 1]] ψz[z, Subscript[Z,i + 1]];  

Note: We use $\mu$ to represent $p$.

Here are the assumptions I make:

BasisAssumY = {Δy > 0, Subscript[Y, i] ∈ Reals,
y ∈ Reals, Δy ∈ Reals, 
Subscript[Y, i - 1] ∈ Reals, 
Subscript[Y, i + 1] ∈ Reals, 
Subscript[Y, i - 1] < Subscript[Y, i] < Subscript[Y, i + 1], 
Subscript[Y, i + 1] - Subscript[Y, i ] == Δy, 
Subscript[Y, i] - Subscript[Y, i - 1 ] == Δy, 
Subscript[Y, i ] - Subscript[Y, i + 1] == -ΔY};

The same assumptions are made in the z direction.

BasisRulesAll = {Subscript[Z, i + 1] - Subscript[Z, 
i ] -> Δz, 
Subscript[Z, i] - Subscript[Z, i - 1 ] -> Δz, 
Subscript[Z, i ] - Subscript[Z, i + 1] -> -Δz, 
Subscript[Y, i + 1] - Subscript[Y, i ] -> Δy, 
Subscript[Y, i] - Subscript[Y, i - 1 ] -> Δy, 
Subscript[Y, i ] - Subscript[Y, i + 1] -> -Δy};

Here is my integrating procedure:

Simplify[Integrate[
Integrate[ψz[z, Subscript[Z, i]]  θ [y, z] μ [y, 
    z], {z, Subscript[Z, i - 1], Subscript[Z, i + 1]}, 
  Assumptions -> BasisAssumZ] ψy[y, Subscript[Y, i]], {y, 
 Subscript[Y, i - 1], Subscript[Y, i + 1]}, 
Assumptions -> 
 BasisAssumY]/(Δz Δy)] //. \BasisRulesAll

This will run for hours without returning a result. In contrast, I have done a similar calculation in 1D (get the notebook here) that takes about 10 seconds.

I did manage to get the integral to evaluate but I had to remove all $y$ terms from the $x$ integral (I grouped together like terms and factored). After doing that, the integral took about 10 seconds to be evaluated, but it was very messy and I cannot be sure that human error was not introduced. It would be simpler if the method above worked in a reasonable amount of time.

Does anyone have any suggestions on how I could improve the performance of the double integration?

Answer 1

After some trial I managed to speed your integration up. The main idea for the optimization is to simplify the symbolic expression as much as possible before throw it in to Integrate. Also, your code can be conciser for many parts of it, I'll mention a few of them in the following part.

This is the expression waiting for simplification:

exp = ψy[y, Y[i]] ψz[z, Z[i]] θ[y, z] μ[y, z]

Notice that I've changed all of your Subscript[a_, b__] in to something like a[b] for Subscript in some case can be quite troublesome so it's better not to use it.

First, since your θ[y, z] and μ[y, z] have the same form and they're in fact just be used as replacing rule, you don't need to define them as 2 separate functions, just use a rule:

functionrule = ϕ_[y, z] -> 
       ϕ[i - 1, j - 1] ψy[y, Y[i - 1]] ψz[z, Z[i - 1]] + 
       ϕ[i - 1, j] ψy[y, Y[i]] ψz[z, Z[i - 1]] + 
       ϕ[i - 1, j + 1] ψy[y, Y[i + 1]] ψz[z, Z[i - 1]] + 
       ϕ[i, j - 1] ψy[y, Y[i - 1]] ψz[z, Z[i]] + 
       ϕ[i, j] ψy[y, Y[i]] ψz[z, Z[i]] +
       ϕ[i, j + 1] ψy[y, Y[i + 1]] ψz[z, Z[i]] +
       ϕ[i + 1, j - 1] ψy[y, Y[i + 1]] ψz[z, Z[i - 1]] +
       ϕ[i + 1, j] ψy[y, Y[i]] ψz[z, Z[i + 1]] +
       ϕ[i + 1, j + 1] ψy[y, Y[i + 1]] ψz[z, Z[i + 1]];

Then I noticed that your BasisAssumY shows that there's a simple relation between X_[i_] and X_[i_+1_], so I think it may help if I introduce these relations as replacing rules rather than assumptions for Integrate:

deltarule = {a_[b_ + c_Integer] :> a[b] + c Symbol["Δ" <> ToLowerCase@ToString@a]}

You notebook for the 1D case and a close look for the definition of your ψz and ψy told me that I can even eliminate the Y_[i] if I shift the domain of Integrate:

shiftrule = {y -> y + Y[i], z -> z + Z[i]};

And finally, after some observation for the result given by the rules above, I think a simplification for those Piecewise function will help, too:

simplifyrule = 
  HoldPattern@Piecewise[a__] :> 
      Simplify[Piecewise[a], Δy > 0 && Δz > 0 && -Δy < y < Δy && -Δz < z < Δz];

OK, plug all this in:

exp = exp /. functionrule /. deltarule /. shiftrule /. simplifyrule

Now, the remaining assumptions are just:

BasisAssum = {Δy > 0, y ∈ Reals, Δz > 0, z ∈ Reals}

Finally, integrate with these assumptions:

(* Notice I've shifted the domain for integration *)
1/(Δz Δy)Integrate[exp, {y, -Δy, Δy}, {z, -Δz, Δz}, 
                   Assumptions -> BasisAssum] // AbsoluteTiming

I use a 2GHz laptop.