Geting error while using quadgk for quadruple integral
显示 更早的评论
I am trying to evaluate quadruple integral using quadgk. It is giving me error even with just one integral. I am pasting the code below
a = 1;
b = 1;
r = 20;
intgrl = quadgk(@(theta1)quadgk(@(theta2)quadgk(@(w)quadgk(@(z)qfunc((a./w.^2)+(b./z.^2)).*(w.^2+z.^2-(4*r^2)<=0).*(w.^2+z.^2-2.*w.*z.*cos(theta1-theta2)>=0).*w.*z.*(1/16*pi^2*r^4), r, 2*r), 0, r), 0,2*pi), 0, 2*pi);
And I am getting the following error.
"Error using + Matrix dimensions must agree."
I feel there is dimension mismatch when evaluating over different variables. Could someone let me know how to evaluate this quadruple integral?
采纳的回答
Mike Hosea
2014-8-11
If you want to integrate a function of 4 variables, f(theta1,theta2,w,z), try iterating it as a single integral of a triple integral like so.
Q1 = @(theta1)integral3(@(theta2,w,z)f(theta1,theta2,w,z),theta2min,theta2max,wmin,wmax,zmin,zmax);
Now you can polish it off by
intgrl = integral(Q1,theta1min,theta1max,'ArrayValued',true);
The 'ArrayValued',true part ensures that integral() will not try to evaluate its integrand function, Q1, at anything other than scalar values of theta1. I'm not saying that this is the fastest way, but the 'ArrayValued' option is the easiest way because otherwise you have to do the work to define your integrand function, f, in such a way that it satisfies the requirements of integral or quadgk, namely that if given arrays of the same size as inputs for all its arguments, it returns an array of the same size. But quadgk returns a scalar, not an array, so you would need to write loops or use arrayfun to define your integrand function, and it can get complicated, or at least complicated-looking.
2 个评论
Sri
2014-8-18
Mike Hosea
2014-9-2
Accuracy-wise, QUADV is not reliable enough for this problem, I think. I don't trust it.
The 'tiled' method of INTEGRAL2 and INTEGRAL3 hates masking functions. In fact masking functions are always a problem numerically, but the 'tiled' method just can't handle them. It has to do with the way it refines in 2D rather than just in 1D. You can switch to the 'iterated' method of INTEGRAL3, though you will probably need to loosen the tolerances substantially for performance reasons. For best results, you should try to reformulate this so that you are integrating qfunc((a./ds.^v)-(a./di.^v))) over the correctly defined domain instead of integrating over a rectangular region and masking it off to be zero over part of the rectangle. The outer boundary of the innermost double integral seems to be the circle ds^2 + di^2 <= (2*r)^2, but the inner boundary is a little more complicated. Bottom line is that defining the region of integration using function handles should work a LOT better if you can sort it all out. It might require splitting the integral up, but it is so expensive to integrate over the discontinuities that the masking functions create, you'll find it to be a bargain.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
