Faster symbolic integration

Question

I have the following three functions $\lambda_i$ and $q_i$:

λ1 = (a x + b y + c)/.First@Solve[{a x1 + b y1 + c == 1, a x2 + b y2 + c == 0, 
     a x3 + b y3 + c == 0}, {a, b, c}];
λ2 = (a x + b y + c)/.First@Solve[{a x1 + b y1 + c == 0, a x2 + b y2 + c == 1, 
     a x3 + b y3 + c == 0}, {a, b, c}];
λ3 = (a x + b y + c)/.First@Solve[{a x1 + b y1 + c == 0, a x2 + b y2 + c == 0, 
     a x3 + b y3 + c == 1}, {a, b, c}];

q1 = 1;
q2 = x - (x1+x2+x3)/3;
q3 = y - (y1+y2+y3)/3;

Then we have 7 functions $\phi_i$ defined as follows

φ1 = λ1 (2 λ1 - 1) + 3 λ1 λ2 λ3;
φ2 = λ2 (2 λ2 - 1) + 3 λ1 λ2 λ3;
φ3 = λ3 (2 λ3 - 1) + 3 λ1 λ2 λ3;
φ4 = 4 λ2 λ3 - 12 λ1 λ2 λ3;
φ5 = 4 λ3 λ1 - 12 λ1 λ2 λ3;
φ6 = 4 λ1 λ2 - 12 λ1 λ2 λ3;
φ7 = 27 λ1 λ2 λ3;

I am interested in integrals in this form $$ \int \limits_{\text{triangle}} \phi_i \: \mathrm{d} S\\ \int \limits_{\text{triangle}} q_i \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} q_i \partial_x \phi_j \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} q_i \partial_y \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \partial_x \phi_j \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \partial_y \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \phi_j \partial_x \phi_k \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \phi_j \partial_y \phi_k \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \partial_x \phi_j \partial_x \phi_k \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \partial_x \phi_j \partial_y \phi_k \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \partial_y \phi_j \partial_y \phi_k \: \mathrm{d} S\\ $$ where "triangle" is defined by its vertices $(x_1, y_1)$, $(x_2, y_2)$, $(x_3, y_3)$

So I tried it with $\phi_i, \phi_j$. It is always proportional to the are of the triangle, so I am interested only in the numeric constant in front of it (integral = something*expression for the area)

Δ = 
  1/2 Abs[x2 y1 - x3 y1 - x1 y2 + x3 y2 + x1 y3 - x2 y3];
resij = {};
For[i = 1, i <= 7, i++,
 For[j = 1, j <= 7, j++,
  expr = ToExpression[ToString[φ] <> ToString[i]]*
     ToExpression[ToString[φ] <> ToString[j]] // Simplify;
  res = 1/Δ Integrate[
       expr, {x, y} ∈ 
        Triangle[{{x1, y1}, {x2, y2}, {x3, y3}}]] // Simplify // 
    Expand;
  Print["(" <> ToString[i] <> "," <> ToString[j] <> ")   " <> 
     ToString[N@res]]
   AppendTo[resij, res];
  ]]

Although it takes a very long time to compute, it gives the desired result:

(1,1)   0.0321429
(1,2)   -0.00674603
(1,3)   -0.00674603
(1,4)   -0.00793651
...

2520 Partition[resij, 7] // MatrixForm

$$ \begin{pmatrix} 81 & -17 &-17 & -20 & 36 & 36 & 27 \\ -17 & 81 & -17 & 36 &-20 & 36 & 27 \\ -17 & -17 & 81 & 36 & 36 & -20 & 27 \\ -20 & 36 & 36 &208 &-16 &-16 & 108 \\ 36 & -20 & 36 & -16 & 208 &-16 & 108 \\ 36 & 36 & -20 & -16 & -16 & 208 & 108 \\ 27 & 27 & 37 & 108 & 108 & 108 & 729 \end{pmatrix} $$

Then I tried it with $\phi_i \phi_j \partial_x \phi_k$, which is always in form $$ \frac{\Delta}{|\Delta|} \left( a y_1 + b y_2 + c y_3 \right) $$ where $\Delta$ is signed triangle area. So I tried:

Δ = 
  1/2 Cross[{x2 - x1, y2 - y1, 0}, {x3 - x1, y3 - y1, 0}][[3]];
resijdxk = {};
For[i = 1, i <= 7, i++,
 For[j = 1, j <= 7, j++,
  For[k = 1, k <= 7, k++,
   expr = (Abs[Δ]/Δ)^-1 ToExpression[
       ToString[φ] <> ToString[i]]*
      ToExpression[ToString[φ] <> ToString[j]]* 
      D[ToExpression[ToString[φ] <> ToString[k]], x] // 
     Simplify;
   res = Integrate[
       expr, {x, y} ∈ 
        Triangle[{{x1, y1}, {x2, y2}, {x3, y3}}]] // Simplify // 
     Expand;
   AppendTo[
    resijdxk, {Coefficient[res, y1], Coefficient[res, y2], 
     Coefficient[res, y3]}];
   Print["(" <> ToString[i] <> "," <> ToString[j] <> "," <> 
     ToString[k] <> ")   " <> "(" <> 
     ToString[N@Coefficient[res, y1]] <> "," <> 
     ToString[N@Coefficient[res, y2]] <> "," <> 
     ToString[N@Coefficient[res, y3]] <> ")"]
   ]]]

(1,1,1)   (0.,0.0309524,-0.0309524)
(1,1,2)   (0.00892857,-0.00297619,-0.00595238)
...

which gets the result, however, one in 40 seconds. It would take 7*7*7*40 seconds = 3 hours and 48 minutes to compute.

What is a faster way to compute this? I have a suspicion that this is just a bunch of spagetti code I wrote and it can be somehow optimalized, but I don't know how.

Answer 1

This his how we you can compute the requestest integrals on the standard triangle {{0,0},{1,0},{0,1}} as follows:

We start by defining the basis functions

λ1 = x;
λ2 = y;
λ3 = 1 - x - y;
q1 = 1;
q2 = x - 1/3;
q3 = y - 1/3;
φ = {λ1 (2 λ1 - 1) + 3 λ1 λ2 λ3, λ2 (2 λ2 - 1) + 3 λ1 λ2 λ3, λ3 (2 λ3 - 1) + 3 λ1 λ2 λ3, 4 λ2 λ3 - 12 λ1 λ2 λ3, 4 λ3 λ1 - 12 λ1 λ2 λ3, 4 λ1 λ2 - 12 λ1 λ2 λ3, 27 λ1 λ2 λ3};
Dφ = D[φ, {{x, y}, 1}];
Φ = {Transpose[{φ}], Dφ, {{q1}, {q2}, {q3}}};

Memoized way to compute integrals of monomials:

ClearAll[int];
mem : int[{a_, b_}] := mem=Integrate[Integrate[x^a y^b, {y, 0, 1 - x}], {x, 0, 1}]

The working horse function; it takes a integer list (expecting to contain1, 2, and 3, only), reads the respective tensors from Φ and computes the integrals of their tensor products.

getIntegrals[idx_?VectorQ] := Module[{a},
  a = 
    Outer[TensorProduct, 
     Sequence @@ Table[Φ[[i]], {i, idx}], 1];
  Developer`ToPackedArray[
   Map[
    Total[
      KeyValueMap[
       {key, value} \[Function] int[key] value,
       Association@CoefficientRules[#, {x, y}]
       ]
      ] &,
    a,
    {ArrayDepth[a]}
    ]
   ]
  ]

Computing all requested integrals on the standard triangle:

results = {
    getIntegrals[{1}],
    getIntegrals[{3, 1}],
    getIntegrals[{3, 2}],
    getIntegrals[{1, 1}],
    getIntegrals[{1, 2}],
    getIntegrals[{1, 1, 2}],
    getIntegrals[{1, 2, 2}]
    }; // AbsoluteTiming // First

1.11174

Here result[[i]] contains the integrals in the i-th row of the OP. For example, the mass matrix (fourth row) can be found as follows:

Flatten[results[[4]], {{1, 3}, {2, 4}}] // TeXForm

$$ \left( \begin{array}{ccccccc} \frac{9}{560} & -\frac{17}{5040} & -\frac{17}{5040} & -\frac{1}{252} & \frac{1}{140} & \frac{1}{140} & \frac{3}{560} \\ -\frac{17}{5040} & \frac{9}{560} & -\frac{17}{5040} & \frac{1}{140} & -\frac{1}{252} & \frac{1}{140} & \frac{3}{560} \\ -\frac{17}{5040} & -\frac{17}{5040} & \frac{9}{560} & \frac{1}{140} & \frac{1}{140} & -\frac{1}{252} & \frac{3}{560} \\ -\frac{1}{252} & \frac{1}{140} & \frac{1}{140} & \frac{13}{315} & -\frac{1}{315} & -\frac{1}{315} & \frac{3}{140} \\ \frac{1}{140} & -\frac{1}{252} & \frac{1}{140} & -\frac{1}{315} & \frac{13}{315} & -\frac{1}{315} & \frac{3}{140} \\ \frac{1}{140} & \frac{1}{140} & -\frac{1}{252} & -\frac{1}{315} & -\frac{1}{315} & \frac{13}{315} & \frac{3}{140} \\ \frac{3}{560} & \frac{3}{560} & \frac{3}{560} & \frac{3}{140} & \frac{3}{140} & \frac{3}{140} & \frac{81}{560} \\ \end{array} \right) $$