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?
θ[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 thatIntegrate
has a syntax for multiple integral e.g.Integrate[Sin[x y], {x, 0, 1}, {y, 0, x}]
? $\endgroup$BasisAssumY
i.e. theΔY
should beΔy
. $\endgroup$