Numerical 2D Integration in MATLAB - Error

Hi All,

I am having issues with using 2d numerical integration in Matlab 2011b.

I am trying to create a function that takes an input compliance matrix (6x6) and outputs an orientational averaged compliance matrix (6x6).

This is done through the use of transformation equations (which are quite complicated).

Ultimately, I take 2 6x6 transformation matrices (full of cos(x), cos(y), sin(x), sin(y)) and multiply them together with the input stiffness matrix, and then attempt to loop through each cell and numerically integrate to get the orientational average.

I have attached my code:

function[S_prime] = orientation_random_2(S)

% 3D Orientational average of compliance matrix as defined by RUC subject

% to periodic BC's

% T1 = Stress Transformation Matrix

% T2 = Strain Transformation Matrix

% S = Compliance Matrix of RUC

% theta = angle from X1 axis

% phi = angle from X3 axis

% Input 6x6 compliance matrix
% Output 6x6 orientational averaged compliance matrix

syms('x','y');

m = cos(y);

n = sin(y);

p = cos(x);

q = sin(x);

pi = 3.14159265359;

S_prime_variable = zeros(6);

S_Prime_Temp = zeros(6);

S_Prime = zeros(6);

T1 = [m^2*p^2,n^2*p^2,q^2,-2*n*p*q,-2*m*p*q,2*m*n*p^2;n^2,m^2,0,0,0,2*m*n;m^2*q^2,-n^2*q^2,p^2,-2*n*p*q,(m+n)*p*q,-2*m*n*q^2;n*m*q,-n*m*q,0,m*p,n*p,(m^2-n^2)*q;m^2*p*q,-n^2*p*q,-q*p,n*(p^2-q^2),m*(p^2-q^2),0;m*n*p,m*n*p,0,-m*q,-n*q,(m^2+n^2)*p];

T2 = [m^2*p^2,n^2*p^2,q^2,-n*p*q,-m*p*q,m*n*p^2;n^2,m^2,0,0,0,m*n;m^2*q^2,-n^2*q^2,p^2,-n*p*q,0.5*(m+n)*p*q,-m*n*q^2;2*n*m*q,-2*n*m*q,0,m*p,n*p,(m^2-n^2)*q;2*m^2*p*q,-2*n^2*p*q,-2*q*p,n*(p^2-q^2),m*(p^2-q^2),0;2*m*n*p,2*m*n*p,0,-m*q,-n*q,(m^2+n^2)*p];

S_prime_variable =  T2\S*T1;

for j=1:6

       for i=1:6

          fun = @(x,y)S_prime_variable(i,j)*sin(x);

          S_Prime(i,j)=quad2d(fun,0,pi,0,2*pi);

      end
  end

and my error:

>> Y = orientation_random_2(X)

Undefined function 'max' for input arguments of type 'sym'.

Error in quad2d/tensor (line 353)

              if any(any(abs(Z1-Z0) > max(ATOL,RTOL*max(abs(Z1),abs(Z0)))))

Error in quad2d (line 164)

[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in orientation_random_2 (line 40)

          S_Prime(i,j)=quad2d(fun,0,pi,0,2*pi);

Any help with this would be greatly appreciated!!

Answer 1

[Walter Roberson]

When you use

fun = @(x,y)S_prime_variable(i,j)*sin(x);

then the value of the arguments "x" and "y" are not substituted into S_prime_variable(i,j). You should use

fun = @(X,Y) double( subs(S_prime_variable(i,j), {x,y}, {X,Y}) ) * sin(X);

You could also consider using matlabFunction()

[Tyler]

Thanks for responding!

When i switched the way the function was defined as you recommended, i got a new error:

Error using * Inner matrix dimensions must agree.

Error in orientation_random_2/@(X,Y)double(subs(S_prime_variable(i,j),{x,y},{X,Y}))*sin(X)

Error in quad2d/tensor (line 346) Z1 = FUN(X(VTSTIDX),Y(VTSTIDX)); NFE = NFE + 1;

Error in quad2d (line 164) [Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in orientation_random_2 (line 37) S_Prime(i,j)=quad2d(fun,0,pi,0,2*pi);

When using the matlabFunction() command, i get this error:

Warning: Reached the maximum number of function evaluations (2000). The result fails the global error test.

> In quad2d at 241

In orientation_random_2 at 39

Warning: Reached the maximum number of function evaluations (2000). The result fails the global error test.

> In quad2d at 241

In orientation_random_2 at 39

Warning: Reached the maximum number of function evaluations (2000). The result fails the global error test.

> In quad2d at 241

In orientation_random_2 at 39

Any suggestions?