# Calculating double integrals from a very big function

I need to compute integrals (for different value of i and j) of a function named se in this post. We need to double differentiate from se related to i and j and then compute double integrals over the domains of x and y. Integrating from se functions are very time consuming and produce errors that makes one doubt whether the results are correct or not. se function is build from chebyshev polynomial double series functions. chebyshev series is defined as below: Where: Uij, Vij, Wij, Xij and Yij are constants for each value of i and j differentiated the se functions according to them.

Ru(x,y), Rv(x,y), Rw(x,y), Rx(x,y) and Ry(x,y) (that are behind the series) are the power of p, q ,r ans s are assumed to be one. It is assumed Nx=Ny and defined by Np in the Mathematica code.

Let start to define the se function. We need some preliminaries. at first, we define some constants of our problem.

 a = 0.0762;
b = 0.1016;
r = 0.762;
b = b/r;
h = 0.0003302;
Em = 68.947*10^9;
ru0 = 2657.27;
nu = 0.33;
G = Em/(2*(1 + nu));

Np = 10;


(Np is the number of series in chebyshev functions))

 Nip = 100;
(*Nip is the number of dividing part of x and y in integration region)*)

Q11 = Em/(-nu^2 + 1);
Q12 = Em*nu/(-nu^2 + 1);
Q22 = Q11;
Q44 = G;
Q55 = G;
Q66 = G;

A11 = Q11*h; A12 = Q12*h; A16 = 0; A22 = Q22*h; A26 = 0; A66 =
Q66*h; A44 = Q44*h; A55 = Q55*h; A45 = 0;

D11 = 1/3*((2*h^3)/8)*Q11;
D12 = Q12*1/3*((2*h^3)/8);
D16 = 0;
D22 = Q22*1/3*((2*h^3)/8);
D26 = 0;
D66 = Q66*1/3*((2*h^3)/8);


I define the constants of Uij, Vij, Wij, Xij and Yij as below and define Dis functions as a combination of them in a matrix with one row:

 DU := Array[U, {Np, Np}];
DV := Array[V, {Np, Np}];
DW := Array[W, {Np, Np}];
DX := Array[X, {Np, Np}];
DY := Array[Y, {Np, Np}];
Dis := Join[ArrayReshape[DU, {1, Np^2}], ArrayReshape[DV, {1, Np^2}],
ArrayReshape[DW, {1, Np^2}], ArrayReshape[DX, {1, Np^2}],
ArrayReshape[DY, {1, Np^2}], 2];


I define the chebyshev functions as below:

  pii[i_] := Cos[(i - 1)*ArcCos[(2*x/a)]];
pjj[j_] := Cos[(j - 1)*ArcCos[(2*y/b)]];

Ru = (1 + 2*x/a)*(1 - 2*x/a)*(1 + 2*y/b)*(1 - 2*y/b);
Rv = (1 + 2*x/a)*(1 - 2*x/a)*(1 + 2*y/b)*(1 - 2*y/b);
Rw = (1 + 2*x/a)*(1 - 2*x/a)*(1 + 2*y/b)*(1 - 2*y/b);
Rx = (1 + 2*x/a)*(1 - 2*x/a)*(1 + 2*y/b)*(1 - 2*y/b);
Ry = (1 + 2*x/a)*(1 - 2*x/a)*(1 + 2*y/b)*(1 - 2*y/b);

up[i_, j_] := Ru*DU[[i, j]]*pii[i]*pjj[j];
vp[i_, j_] := Rv*DV[[i, j]]*pii[i]*pjj[j];
wp[i_, j_] := Rw*DW[[i, j]]*pii[i]*pjj[j];
fixp[i_, j_] := Rx*DX[[i, j]]*pii[i]*pjj[j];
fiyp[i_, j_] := Ry*DY[[i, j]]*pii[i]*pjj[j];

u = Sum[up[z1, z2], {z1, Np}, {z2, Np}];
v = Sum[vp[z1, z2], {z1, Np}, {z2, Np}];
w = Sum[wp[z1, z2], {z1, Np}, {z2, Np}];
fix = Sum[fixp[z1, z2], {z1, Np}, {z2, Np}];
fiy = Sum[fiyp[z1, z2], {z1, Np}, {z2, Np}];


After some differentiations and algebraic computations as below:

 exx0 = D[u, x];
eyy0 = D[v, y] + w/r;
gamaxy0 = D[v, x] + D[u, y];
Xxx = D[fix, x];
Xyy = D[fiy, y];
Xxy = D[fix, y] + D[fiy, x];
gamaxz = D[w, x] + fix;
gamayz = D[w, y] + fiy;

Nx = A11*exx0 + A12*eyy0 + A16*gamaxy0;
Ny = A12*exx0 + A22*eyy0 + A26*gamaxy0;
Nxy = A16*exx0 + A26*eyy0 + A66*gamaxy0;
Mx = D11*Xxx + D12*Xyy + D16*Xxy;
My = D12*Xxx + D22*Xyy + D26*Xxy;
Mxy = D16*Xxx + D26*Xyy + D66*Xxy;
Qx = 5/6*(A45*gamayz + A55*gamaxz);
Qy = 5/6*(A44*gamayz + A45*gamaxz);


We define se function as:

 se = 0.5*r*(Nx*exx0 + Ny*eyy0 + Nxy*gamaxy0 + Mx*Xxx + My*Xyy +
Mxy*Xxy + Qy*gamayz + Qx*gamaxz);


So, it the time to define integrals. we should to take the double integrals from the second derivative of the se function with respect to the constant coefficients in Dis matrix. I defined this double integrals as follow according to trapezoidal rule of integration:

 K1[i0_, j0_] := Module[{i = i0, j = j0},
ix1 = Table[{x, D[se, Dis[[1, i]], Dis[[1, j]]]}, {x, -a/2, a/2,a/Nip}];
ix2 = 1/2*Total[Differences[ix1[[All, 1]]] ListCorrelate[{1, 1},ix1[[All, 2]]]];

iy1 = Table[{y, ix2}, {y, -b/2, b/2, b/Nip}];
iy2 = 1/2*Total[Differences[iy1[[All, 1]]] ListCorrelate[{1, 1},iy1[[All, 2]]]]];


At last, we need to define the matrix of L1 with i product j dimensions (because we need to differentiate from se related to each i and j, so a i product j matrix is produced in this procedure) that i and j are 5*Np^2, as below;

 proc1[n_, f_] := Block[{res}, res = ConstantArray[0, {n, n}];
Do[res[[i, j]] = res[[j, i]] = K1[i, j], {i,
n}, {j, i}]; res];

L1 = proc1[5*Np^2, Lm];// Timing


L1 is a matrix with 5Np^2 product 5Np^2 dimensions that we need as results. For Np is 1, 21 second is elapsed for program, and For Np=2, 95 seconds is elapsed (my computer is old, maybe in another pc, fewer time will be displayed). I need to consider Np 14 or more, but with this procedure, for Np= 4 or 5, The process takes several hours and because of errors, I am not sure the results are true? I will be appreciated to guide me to fastening the speed of computations and Overcoming the errors.

Thank you

• You can probably gain some speed by learning when to use immediate assignments = and when to use delayed assignments :=. You are delaying many of the assignments unnecessarily, which means that a lot of your code is executed repeatedly where it could be done once only. wolfram.com/language/elementary-introduction/2nd-ed/… Mar 8 at 10:40
• thank you Roman for your answer. I tried your way and replaced some := with = that I think help to improve the speed of code, but the results faced with troubles. I tested different situations in replacing but couldn't get the results. Could you please see my code and tell which one of := can be changeable with = ? Mar 8 at 17:21