Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
Do you mean that the θ[y, z] and μ[y, z] appearing at the end of the code take the form of your ϕ[y, z] by saying "The p term takes the same form but is represented by μ"? Also, are you aware that Integrate has a syntax for multiple integral e.g. Integrate[Sin[x y], {x, 0, 1}, {y, 0, x}]? –  xzczd Feb 27 at 13:02
    
I added the function for $\mu$ in. I accidentally put the $\phi$ function instead of our $\theta$ function. I did know that you could integrate using a single command but it does not speed things up. –  rioneye Feb 27 at 15:47
    
Based on the context of your question, I think there's a typo in the last element of your BasisAssumY i.e. the ΔY should be Δy. –  xzczd Feb 28 at 4:12
add comment

1 Answer 1

up vote 1 down vote accepted

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

enter image description here

I use a 2GHz laptop.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.