Improve the speed of Gaussian quadrature integration

Question

I am using Gaussian quadrature method to do numerical integration for a function Pderivative, which is related to the demagnetization tensor of elastic-electric-magnetic coupling field. Pderivative is a 17*17 matrix function built by the symbolic partial derivative of a function stored in a file ptnsor2.m. The code is shown below.

Needs["NumericalDifferentialEquationAnalysis`"]
P1 = << "ptensor2.m";
P2[c11_, c12_, c13_, c33_, c44_, e31_, e33_, e15_, q31_, q33_, q15_, 
k11_, k33_, a11_, a33_, u11_, u33_, t_] = {P1[[1, 1]], P1[[1, 2]], 
P1[[1, 3]], P1[[3, 3]], P1[[4, 4]], P1[[1, 9]], P1[[3, 9]], 
P1[[5, 7]], P1[[1, 12]], P1[[3, 12]], P1[[5, 10]], P1[[7, 7]], 
P1[[9, 9]], P1[[7, 10]], P1[[9, 12]], P1[[10, 10]], P1[[12, 12]]};
Pderivative[c11_, c12_, c13_, c33_, c44_, e31_, e33_, e15_, q31_, 
q33_, q15_, k11_, k33_, a11_, a33_, u11_, u33_, t_] = 
D[P2[c11, c12, c13, c33, c44, e31, e33, e15, q31, q33, q15, k11, 
k33, a11, a33, u11, u33, 
t], {{c11, c12, c13, c33, c44, e31, e33, e15, q31, q33, q15, k11, 
 k33, a11, a33, u11, u33}}];
wt = GaussianQuadratureWeights[20, -1, 1];
Ptensor[A_] := 
Module[{c11, c12, c13, c33, c44, e31, e33, e15, q31, q33, q15, k11, 
k33, a11, a33, u11, u33},
c11 = A[[1, 1]]; c12 = A[[1, 2]]; c13 = A[[1, 3]]; c33 = A[[3, 3]]; 
c44 = A[[4, 4]];
e31 = A[[1, 9]]; e33 = A[[3, 9]]; e15 = A[[5, 7]];
q31 = A[[1, 12]]; q33 = A[[3, 12]]; q15 = A[[5, 10]];
k11 = -A[[7, 7]]; k33 = -A[[9, 9]];
a11 = -A[[7, 10]]; a33 = -A[[9, 12]];
u11 = -A[[10, 10]]; u33 = -A[[12, 12]];
Table[Sum[
Pderivative[c11, c12, c13, c33, c44, e31, e33, e15, q31, q33, q15,
    k11, k33, a11, a33, u11, u33, wt[[k, 1]]][[i, j]]*
 wt[[k, 2]], {k, 1, 20}], {i, 1, 17}, {j, 1, 17}]]

The file ptensor2.m can be accessed via the following link: https://drive.google.com/file/d/0By4JUh1c-zmbMkRvTTNvenp5OFk/view?usp=sharing

I am using Thinkpad T420 with Intel i7-2640M CPU, and mathematica V10 on Window 7 system. I try to run this integration for a specific input matrix, and it turns out to take a long time.

A0 = {{1.6603956637004407`*^11, 7.702892284413217`*^10, 
7.80285626281322`*^10, 0.`, 0.`, 0.`, 0.`, 
0.`, -4.397908045934428`, 0.`, 0.`, 
0.023595669736661423`}, {7.702892284413217`*^10, 
1.6603956637004407`*^11, 7.80285626281322`*^10, 0.`, 0.`, 0.`, 
0.`, 0.`, -4.397908045934428`, 0.`, 0.`, 
0.023595669736661478`}, {7.80285626281322`*^10, 
7.80285626281322`*^10, 1.620486356879899`*^11, 0.`, 0.`, 0.`, 0.`,
 0.`, 18.585111713980364`, 0.`, 0.`, 0.03440290158172471`}, {0.`, 
0.`, 0.`, 4.300379639073157`*^10, 0.`, 0.`, 0.`, 
11.591885622145135`, 0.`, 0.`, 0.008028677439673062`, 0.`}, {0.`, 
0.`, 0.`, 0.`, 4.300379639073157`*^10, 0.`, 11.591885622145135`, 
0.`, 0.`, 0.008028677439673062`, 0.`, 0.`}, {0.`, 0.`, 0.`, 0.`, 
0.`, 4.450532176295594`*^10, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`}, {0.`, 
0.`, 0.`, 0.`, 11.591885622145135`, 0.`, -1.1192610752980372`*^-8,
 0.`, 0.`, 6.985516113011081`*^-13, 0.`, 0.`}, {0.`, 0.`, 0.`, 
11.591885622145135`, 0.`, 0.`, 0.`, -1.1192610752980372`*^-8, 0.`,
 0.`, 6.985516113011081`*^-13, 
0.`}, {-4.397908045934428`, -4.397908045934428`, 
18.585111713980364`, 0.`, 0.`, 0.`, 0.`, 
0.`, -1.2590850087609594`*^-8, 0.`, 
0.`, -3.747136660265141`*^-13}, {0.`, 0.`, 0.`, 0.`, 
0.008028677439672059`, 0.`, 6.985516112928487`*^-13, 0.`, 
0.`, -5.008506296673207`*^-6, 0.`, 0.`}, {0.`, 0.`, 0.`, 
0.008028677439672059`, 0.`, 0.`, 0.`, 6.985516112928487`*^-13, 
0.`, 0.`, -5.008506296673207`*^-6, 0.`}, {0.023595669736661138`, 
0.023595669736661138`, 0.034402901581718125`, 0.`, 0.`, 0.`, 0.`, 
0.`, -3.7471366602224796`*^-13, 0.`, 
0.`, -0.000010010069026456233`}};
Timing[result = Ptensor[A0];]

The Timing result is 929.781560, which is too long for me. I plan to do about 4000 to 5000 iterations with the above calculation, so I hope the run time could be reduced to less than 10s. Is there any way I can achieve this goal, like optimizing codes, or running on a cluster (if running on a cluster, please specify if the codes needs be modified first)? Thank you.

Well one solution will be trying to compile Pderivative while using Internal\Bag` to store matrices. Check here: http://mathematica.stackexchange.com/questions/845/internalbag-inside-compile However, given the complexity of your code, I am afraid this will be a tedious task. Can you give more details on how you actually obtained ptensor2? Maybe we can optimize from there as such calculation will definitely take long time .. — Bichoy, Jun 15 '15 at 06:00
Getting a full symbolic expression then substituting numbers is not always the best route especially if the expressions becomes so complicated. Dealing with the problem as a numerical one from the beginning might be a more optimal solution. Think about inverting matrices as one example where direct application of numerical method is super fast compared to evaluating analytical expressions. — Bichoy, Jun 15 '15 at 06:03
@Bichoy I finally solve this problem by using Matlink to make this part of code run in Matlab and then return the results to Mathematica. It turns out Matlab is pretty fast, only taking less than 2s. I don't know why it takes Mathematica so long... But anyway thanks for the comments and discussions. — riorita, Jun 17 '15 at 18:21
I wouldn't be surprised to learn that MATLAB is faster when dealing with matrices ... However, I still believe we could do a much better job if from the beginning, we avoid generating this large symbolic expressions ... — Bichoy, Jun 19 '15 at 02:26

Improve the speed of Gaussian quadrature integration

0 Answers0