# Faster symbolic integration

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.

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];
DeveloperToPackedArray[
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)$$

• The first solution looks fine (giving the desired result, up to a constant). However, I don't get the second solution. I'm looking for coefficients that appear in front of the x1, x2, x3 or y1, y2, y3 when I divide that integral by the sign of signed integral area. Your results just give one number per integral. Other than that, I think that Outer[Times, phi, phi, phi] should have the last one D[phi, x] or D[phi, y], but that's just detail. I'm more interested in those constants, rather than the pure value of the integral over the basic simplex. – user16320 Jan 9 '18 at 17:32