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);de here

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 M1 and K1 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 these double integrals as follow according to trapezoidal rule of integration:

M1 := Compile[{i, j},
Block[{x, y, Dis, t, Nip},
ix1 = 
 Table[{x, D[t, 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]]]]; iy2]];


K1 := Compile[{i, j},
Block[{x, y, Dis, se, Nip},
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]]]]; iy2]];

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_] := 
    proc1[n, f] = Block[{res}, res = ConstantArray[0, {n, n}];
      Do[res[[i, j]] = res[[j, i]] = K1[i, j] - omega^2*M1[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 15 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. **

All the above codes (in a notebook file) is in the attached file

Thank you

Attachments:

perfect panel1.nb