Area integral over a triangle in MATLAB -- is numerical integration possible?
in [Matlab]
|
From: Roger Stafford on 22 Mar 2010 11:19 "Vivek Saxena" <maverick280857(a)yahoo.com> wrote in message <ho6g8s$os1$1(a)fred.mathworks.com>... > ......... > Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. ....... ------------ One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)? If so, numerical integration is a correspondingly difficult process and dblquad, or possibly even quad2d, may very well end up confused and erroneous unless you place a very tight tolerance setting in them. This is true whether or not you have performed a transformation of the variables, and it is also true even for single integration. Quadrature routines can be temperamental beasts when faced with rapidly fluctuating integrands. You would do well to read the description of quad2d at http://www.mathworks.com/access/helpdesk/help/techdoc/ref/quad2d.html to get an understanding of these difficulties even if you are only using dblquad. Roger Stafford
From: Bruno Luong on 22 Mar 2010 11:31 "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <ho81p7$1nv$1(a)fred.mathworks.com>... > One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)? Roger, I can answer this question. Usually not: when the FEM is correctly setup then the mesh size must be small compare to wave number, thus exp(-i*k_0*x) would not have many cycles on the triangle element. If it's not the case then there is no use to solve FEM because the discretized solution will have huge error. Bruno
From: Vivek Saxena on 22 Mar 2010 12:01 "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <ho82fo$eg4$1(a)fred.mathworks.com>... > "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <ho81p7$1nv$1(a)fred.mathworks.com>... > > > One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)? > > Roger, I can answer this question. Usually not: when the FEM is correctly setup then the mesh size must be small compare to wave number, thus exp(-i*k_0*x) would not have many cycles on the triangle element. If it's not the case then there is no use to solve FEM because the discretized solution will have huge error. > > Bruno Hi Roger and Bruno. Thanks for your replies and help so far. I'm using a unit wavelength, and k_0 = 2*pi as a result. Also, the mesh size has been kept small enough -- actually this is a textbook problem but I have to solve it with a more complicated boundary condition than is given in the standard books. So these particular double integrals arise in the second order absorbing boundary condition, but not in the formulation with the first order absorbing boundary condition (for which the code runs perfectly). In fact, theory tells me that the second order ABC should work fine with a coarser discretization as compared to the first order ABC (if things are going fine). But I am being conservative and keeping the mesh size just as fine as it was earlier. Roger, I think it is likely that the program is plagued by some numerical errors, as you have pointed out. (I have gone through the code and ensured that there are no logical or syntactic errors.) I must admit that I do not know these techniques very well theoretically, but I will try and play with the options. If you have more suggestions on how the tolerances could be set in view of the kind of function that the integrand is, please do advise me.
From: Roger Stafford on 23 Mar 2010 13:03 "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <ho75n5$b42$1(a)fred.mathworks.com>... > One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of > - order 0: the integral result is P(x1)*|T| > - order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3 > etc... > > x1, x2, x3 are three corners, and |T| is the area of triangle > > Once the implementation gives correct result, you can then replace P(x) by your function. > > Bruno ------------ Hi Bruno. Your "etc..." intrigued me, so I worked out what the "order 2" test would be. I wasn't aware of this identity before. What would an order 3 test be? The double integral of any quadratic polynomial, P(x,y) = A*x^2+B*x*y+C*y^2+D*x+E*y+F , over a triangle is equal to the triangle's area, multiplied by one quarter of the polynomial's average value at the three vertices plus three quarters of its value at the triangle's centroid (located at the averages of the three respective coordinates.) Roger Stafford
From: Bruno Luong on 23 Mar 2010 14:11
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hoas8r$jbd$1(a)fred.mathworks.com>... > "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <ho75n5$b42$1(a)fred.mathworks.com>... > > One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of > > - order 0: the integral result is P(x1)*|T| > > - order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3 > > etc... > > > > x1, x2, x3 are three corners, and |T| is the area of triangle > > > > Once the implementation gives correct result, you can then replace P(x) by your function. > > > > Bruno > ------------ > Hi Bruno. Your "etc..." intrigued me, so I worked out what the "order 2" test would be. I wasn't aware of this identity before. What would an order 3 test be? > > The double integral of any quadratic polynomial, > > P(x,y) = A*x^2+B*x*y+C*y^2+D*x+E*y+F , > > over a triangle is equal to the triangle's area, multiplied by one quarter of the polynomial's average value at the three vertices plus three quarters of its value at the triangle's centroid (located at the averages of the three respective coordinates.) > > Roger Stafford Hi Roger, The order 3 polynomial, the formula is I = 9/20*T*[P((x1+x2+x3)/3)] + 2/15*T*[P((x1+x2)/2)+P((x2+x3)/2)+P((x3+x1)/2)] + 1/20*T*[P(x1) + P(x2) + P(x3)]. Those kind of Gauss's formula are derived mostly for FEM (I copy it from a book). Bruno |