How to Use Trapezoid Method with 3 Loops in MATLAB for Numeric Integration

[quin]

dear users
I want to do a triple integration over x ,y ,z in MATLAB but because my expression under the integral is so complicated I want to change integrals to sum by Trapezoid method using a loop in MATLAB but for doing that I have 3 integrals and so 3 for loops for summing over x , y , z
can you please tell me Trapezoid method using 3 loops in matlab?

 

[kreil]

For a discrete set of points,

\int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)= \frac{b-a}{2N}(f(x_1) + 2f(x_2) + 2f(x_3) + \ldots + 2f(x_N) + f(x_{N+1}))

I wrote a simple MATLAB program a few years ago that implements this for one variable. It should be straightforward to extend this to two or more variables- the process is the same. I attached the m-file.

 

[quin]

For a discrete set of points,

\int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)= \frac{b-a}{2N}(f(x_1) + 2f(x_2) + 2f(x_3) + \ldots + 2f(x_N) + f(x_{N+1}))

I wrote a simple MATLAB program a few years ago that implements this for one variable. It should be straightforward to extend this to two or more variables- the process is the same. I attached the m-file.

thanks a lot but sorry it has some errors that I don't know how to solve them
would you mind cheking your code?
and also what is the form of funhandle that I must write in function? forexample if I write :myTrap(x,[0,1],100) it takes errors.
thank you

 

[kreil]

It takes function handles as inputs. For example, for f(x) = x,

fcn = @(x) x;
myTRAP(fcn,[0 1], 100)

ans =

    0.5000

Or for sin(x)cos(x),

fcn = @(x) sin(x).*cos(x);
myTRAP(fcn,[0 1], 100)

ans =

    0.3540

 

[quin]

It takes function handles as inputs. For example, for f(x) = x,

fcn = @(x) x;
myTRAP(fcn,[0 1], 100)

ans =

    0.5000

Or for sin(x)cos(x),

fcn = @(x) sin(x).*cos(x);
myTRAP(fcn,[0 1], 100)

ans =

    0.3540

theanks a lot but the thing that make my problem hard is that I don't have an explicit function in terms of x,y,z. My function will be made through those sums of integral.
see, I have one matrix that I must do set of algabric operation on its eigenvalues and after finding an expression in terms of x,y,z I must integrate that expression(triplequad)

But finding its eigenvalues is impossible because of complexibility
So I must write a code that at first turn "x,y,z" symbols to 3 "for loops"( so that my algabric operation become numeric) and after that sum the result of those 3 loops in away that make trapezoid algorithm
I hope that I could tell my purpose because I really need help for writing my code please !

 

[kreil]

Ok, I'll see what I can do. The first step has to do with the matrix. Do you want to do all of this for one specific matrix? Perhaps for a specific class of matrices (all 3x3x3s)?

Depending on the conditions for the input matrix, we can probably use the EIG function to find the eigenvalues and eigenvectors.

 

[quin]

No my matrix is defferent in size forexample 12*12 and…
but I now just write the 4*4 matrix that is simple and want to test code
I wrote a code and want to tell my purpose easier
This is all the thing that I want to dos=[ ] 4*4 matrix in terms of x,y,z
s1 , s2,s3,s4 are its eigenvalues
I must do this integral and renges of all of 3 integrals are from pi to -pi:
∫∫∫ (1/(m+S1)) + (1/(m+S2)) + (1/(m+S3)) +(1/(m+S4)) dx dy dz

but because for my main matrix(not what I wrote here) I cannot find eigenvalues I must change it to numerical matrix. so I change integral to 3 for loops ( in my code) and for each x,y,z I find eig and sum over all of possible sentences.
but my answer is not so accurate so I must change the way of integration from rectangle to trapezoid.
thanks a lot for spending time on my problem.

 

[kreil]

You can use symbolic math to do this.

syms A x y z
A = 2 * [1 cos((y+z)/4) cos((x+z)/4) cos((x+y)/4); cos((y+z)/4) 
1 cos((x-y)/4) cos((x-z)/4); cos((x+z)/4) cos((x-y)/4) 1 cos((y-z)/4); 
cos((x+y)/4) cos((x-z)/4) cos((y-z)/4) 1;];

S = eig(A)

S =
 
 2 - (12*((64*cos(x/4)^2*cos(y/4)^2*cos(z/4)^2 + 64*sin(x/4)^2*sin(y/4)^2*sin(z/4)^2)^2/2 - (8*cos(x/4)^2*cos(y/4)^2 + 8*cos(x/4)^2*cos(z/4)^2 + 8*cos(y/4)^2*cos(z/4)^2 + 8*sin(x/4)^2*sin(y/4)^2 + 8*sin(x/4)^2*sin(z/4)^2 + 8*sin(y/4)^2*sin(z/4)^2)^3/27 + (3^(1/2)*(27*(64*cos(x/4)^2*cos(y/4)^2*cos(z/4)^2 + 64*sin(x/4)^2*sin(y/4)^2*sin(z/4)^2)^4 - 16*(8*cos(x/4)^2*cos(y/4)^2 + 8*cos(x/4)^2*cos(z/4)^2 + 8*cos(y/4)^2*cos(z/4)^2 + 8*sin(x/4)^2*sin(y/4)^2 + 8*sin(x/4)^2*sin(z/4)^2 + 8*sin(y/4)^2*sin(z/4)^2)^4*(16*cos(x/4)^4*cos(y/4)^4 + 16*cos(x/4)^4*cos(z/4)^4 + 16*cos(y/4)^4*cos(z/4)^4 + 16*sin(x/4)^4*sin(y/4)^4 + 16*sin(x/4)^4*sin(z/4)^4 + 16*sin(y/4)^4*sin(z/4)^4 - 32*cos(x/4)^2*cos(y/4)^2*cos(z/4)^4 - 32*cos(x/4)^2*cos(y/4)^4*cos(z/4)^2 - 32*cos(x/4)^4*cos(y/4)^2*cos(z/4)^2 - 32*sin(x/4)^2*sin(y/4)^2*sin(z/4)^4 - 32*sin(x/4)^2*sin(y/4)^4*sin(z/4)^2 - 32*sin(x/4)^4*sin(y/4)^2*sin(z/4)^2 - 32*cos(x/4)^2*cos(y/4)^2*sin(x/4)^2*sin(y/4)^2 + 32*cos(x/4)^2*cos(y/4)^2*sin(x/4)^2*sin(z/4)^2 + 32*cos(x/4)^2*cos(z/4)^2*sin(x/4)^2*sin(y/4)^2 + 32*cos(x/4)^2*cos(y/4)^2*sin(y/4)^2*sin(z/4)^2 - 32*cos(x/4)^2*cos(z/4)^2*sin(x/4)^2*sin(z/4)^2 + 32*cos(y/4)^2*cos(z/4)^2*sin(x/4)^2*sin(y/4)^2 + 32*cos(x/4)^2*cos(z/4)^2*sin(y/4)^2*sin(z/4)^2 + 32*cos(y/4)^2*cos(z/4)^2*sin(x/4)^2*sin(z/4)^2 - 32*cos(y/4)^2*cos(z/4)^2*sin(y/4)^2*sin(z/4)^2) - 4*... Output truncated.  Text exceeds maximum line length of 25,000 characters for Command Window display.

S now contains the symbolic answer for the eigenvalues. You can plug in actual values for x, y, z by using subs:

r = subs(S,[x y z],[1 1 1]);
eval(r)

ans =

  -0.0000 - 0.0000i
  -0.0000 + 0.0000i
   0.3611 + 0.0000i
   7.6389 + 0.0000i

 

[quin]

my main matrix that is 12*12 is complex for MATLAB to find its eigenvalues.in above code I just wrote this 4*4 as an example for explaining my purpose.
matlab cannot find its eigenvalues in terms of x,y,z

 

[kreil]

I edited the code a bit and am able to get an answer with a reasonable amount of time now. I used some random values for some undefined variables from your scipt, but you can edit these to the values relevant to you. Here is what I have:

tic;
de = 0.2;
d=1;
F=0;
m=1;
for z=-pi:de:pi
  for x=-pi:de:pi
    for y=-pi:de:pi
j=[(2*d)/3  (2*d)/3    (2*d)/3      -(2*cos((x +z)/4))   0    0        -(2*cos((y +z)/4))  0    0      -(2*cos((x +y)/4))  0   0;
  (2*d)/3    (2*d)/3  (2*d)/3       0   -(2*cos((x +z)/4))   0        0   -(2*cos((y +z)/4))   0     0   -(2*cos((x +y)/4))   0;
  (2*d)/3    (2*d)/3    (2*d)/3     0    0     -(2*cos((x +z)/4))     0    0    -(2*cos((y +z)/4))   0    0     -(2*cos((x +y)/4));
  -(2*cos((x +z)/4))   0    0         (2*d)/3  -(2*d)/3    -(2*d)/3   -(2*cos((x -y)/4))  0   0      -(2*cos((y-z)/4))  0    0;
  0    -(2*cos((x +z)/4))   0        -(2*d)/3    (2*d)/3  (2*d)/3     0   -(2*cos((x -y)/4))  0      0   -(2*cos((y-z)/4))   0;
  0    0     -(2*cos((x +z)/4))      -(2*d)/3    (2*d)/3    (2*d)/3   0    0    -(2*cos((x -y)/4))   0    0     -(2*cos((y-z)/4));
 -(2*cos((y +z)/4))  0   0      -(2*cos((x -y)/4))  0    0        (2*d)/3  -(2*d)/3    (2*d)/3      -(2*cos((x-z)/4))  0    0;
  0   -(2*cos((y +z)/4))  0      0   -(2*cos((x -y)/4))   0        -(2*d)/3    (2*d)/3  -(2*d)/3      0   -(2*cos((x-z)/4))   0;
  0    0    -(2*cos((y +z)/4))   0    0     -(2*cos((x -y)/4))     (2*d)/3    -(2*d)/3    (2*d)/3     0    0    -(2*cos((x-z)/4));
  -(2*cos((x +y)/4))   0   0         -(2*cos((y-z)/4))  0    0      -(2*cos((x-z)/4))  0   0      (2*d)/3  (2*d)/3    -(2*d)/3;
  0   -(2*cos((x +y)/4))   0        0   -(2*cos((y-z)/4))    0     0   -(2*cos((x-z)/4))   0      +(2*d)/3   (2*d)/3  -(2*d)/3;
  0    0     -(2*cos((x +y)/4))     0    0    -(2*cos((y-z)/4))    0    0    -(2*cos((x-z)/4))   -(2*d)/3   -(2*d)/3    (2*d)/3;];
 
S = eig(j,'nobalance');
F = F + 1/(m+S(1))+1/(m+S(2))+1/(m+S(3))+1/(m+S(4))+1/(m+S(5))+1/(m+S(6))+1/(m+S(7))+1/(m+S(8))+1/(m+S(9))+1/(m+S(10))+1/(m+S(11))+1/(m+S(12));
    end
  end
end
Q = F.*(de)^3
toc

The output is

Q =

  551.7377

Elapsed time is 1.340846 seconds.

Note that I haven't double checked any of your formulas here, I just tried to get the program to run (since before it got hung up on eig)

 

[quin]

thanks for your guide
But it seems that you didnt understand what I meant.Let me tell you in another way
forget about that special 12*12 matrix or even 4*4 matrix
we have one matrix in terms of x,y,z, we want its eigenvalues and after finding them( which are still in terms of x,y,z) we want t integrate over a sentence that include those eigenvalues

But because my matrix is complex , Matlab can not find its eigenvalues in symbolic way( I mean if we give x,y,z as an symbol , MATLAB cannot find eigenvalue)
so we change the steps of process:

(instead of first, finding eigenvalues and after that integrate from them) now we first ,change integral to numeric summation and so by this trick , matrix will be numeric and x,y,z isnot symbol any way.and in second step We find eigenvalues

But during change integral to summation we can use different ways like rectangle or trapezoid method.
in the code that I sent you before, I used from rectangle summation(or Riemaan method) that sum each F with last F and after that multiple in small range of sums (de)
But now I want to use trapezois methode.
I hope that I could express my purpose completely.
thank you