Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Trapezoid method using a loop

  1. Mar 3, 2013 #1
    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?
     
  2. jcsd
  3. Mar 3, 2013 #2

    kreil

    User Avatar
    Gold Member

    For a discrete set of points,

    [tex] \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}))[/tex]

    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.
     

    Attached Files:

  4. Mar 4, 2013 #3
    thanks alot but sorry it has some errors that I dont 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
     
  5. Mar 4, 2013 #4

    kreil

    User Avatar
    Gold Member

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

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

    ans =

        0.5000
     
    Or for sin(x)cos(x),

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

    ans =

        0.3540
     
     
  6. Mar 4, 2013 #5
    theanks alot but the thing that make my problem hard is that I dont 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 !
     
  7. Mar 4, 2013 #6

    kreil

    User Avatar
    Gold Member

    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.
     
  8. Mar 4, 2013 #7
    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 do


    s=[ ] 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 alot for spending time on my problem.
     

    Attached Files:

    • trap.m
      trap.m
      File size:
      830 bytes
      Views:
      51
  9. Mar 4, 2013 #8

    kreil

    User Avatar
    Gold Member

    You can use symbolic math to do this.

    Code (Text):

    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:

    Code (Text):

    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
     
     
  10. Mar 4, 2013 #9
    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
     

    Attached Files:

  11. Mar 5, 2013 #10

    kreil

    User Avatar
    Gold Member

    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:

    Code (Text):

    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
    Code (Text):

    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)
     
    Last edited: Mar 5, 2013
  12. Mar 5, 2013 #11
    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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook