Matlab integration error: Integrand output size does not match the input size

28 views (last 30 days)
I am attempting a 4th order nested integration in Matlab. The function is as follows:
second_integral = @(r_dash, phi_dash, r, phi) ( (exp(-1i * k * r_dash * cos(phi_dash) * sin(theta)) .* exp(-1i * k * sqrt(r.^2 + r_dash.^2 - 2*( r * r_dash .* cos (phi_dash - phi))))) / (2 .* pi .* sqrt(r.^2 + r_dash.^2 - 2.*( r .* r_dash .* cos (phi_dash - phi))) ) );
With the constants in this particular example being:
a_1 = 7e-5;
a_2 = 5e-5;
k = 1;
theta = 0; % can also be 0.5, 1, 1.5
The integration is being done as
integrated_second_integral =integral(@(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi)),0,a_1,'ArrayValued',true);
When the integration command is run I get the following errors:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.725290e-19.
> In @(r_dash,phi_dash,r,phi)((exp(-1i*k*r_dash*cos(phi_dash)*sin(theta)).*exp(-1i*k*sqrt(r.^2+r_dash.^2-2*(r*r_dash.*cos(phi_dash-phi)))))/(2.*pi.*sqrt(r.^2+r_dash.^2-2.*(r.*r_dash.*cos(phi_dash-phi)))))
In @(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi)
In integral3>@(y,z)FUN(x(1)*ones(size(z)),y,z) at 139
In funfun/private/integral2Calc>integral2t/tensor at 229
In funfun/private/integral2Calc>integral2t at 56
In funfun/private/integral2Calc at 10
In integral3>innerintegral at 138
In funfun/private/integralCalc>iterateScalarValued at 314
In funfun/private/integralCalc>vadapt at 133
In funfun/private/integralCalc at 76
In integral3 at 122
In @(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi))
In funfun/private/integralCalc>iterateArrayValued at 157
In funfun/private/integralCalc>vadapt at 131
In funfun/private/integralCalc at 76
In integral at 89
In this_matlab_script at 61
Error using integral2Calc>integral2t/tensor (line 242)
Integrand output size does not match the input size.
Error in integral2Calc>integral2t (line 56)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);
Error in integral2Calc (line 10)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);
Error in integral3/innerintegral (line 138)
Q1 = integral2Calc( ...
Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);
Error in integralCalc/vadapt (line 133)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);
Error in integralCalc (line 76)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);
Error in integral3 (line 122)
Q = integralCalc(@innerintegral,xmin,xmax,integralOptions);
Error in
@(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi))
Error in integralCalc/iterateArrayValued (line 157)
fxj = FUN(t(1)).*w(1);
Error in integralCalc/vadapt (line 131)
[q,errbnd] = iterateArrayValued(u,tinterval,pathlen);
Error in integralCalc (line 76)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);
Error in integral (line 89)
Q = integralCalc(fun,a,b,opstruct);
Error in this_matlab_script (line 61)
integrated_second_integral
=integral(@(r_dash)integral3(@(phi_dash,r,phi)second_integral(r_dash,phi_dash,r,phi),0,(2*pi),0,a_2,0,(2*pi)),0,a_1,'ArrayValued',true);
The source of this error is a complete mystery to me. Could someone kindly shed some light on what may be happening here and what I can do to remedy it?

Answers (1)

Mike Hosea
Mike Hosea on 27 Sep 2014
When you see that error, it means that you used a / instead of ./ . The integrator passed in an array for "vectorized" element-wise evaluation, and instead of just processing it elementwise, your integrand ended trying to solve a linear system with it. This problem is a little bit on the edge of what the routines can solve. I had some problems using the default tolerances, so I had to loosen them a bit and then see how far I could tighten them back up.. After fixing the integrand to be vectorized in all variables, I got this
a_1 = 7e-5;
a_2 = 5e-5;
k = 1;
theta = 0;%can also be 0.5,1,1.5
second_integral=@(r_dash,phi_dash,r,phi)((exp(-1i*k*r_dash.*cos(phi_dash)*sin(theta)).*exp(-1i*k*sqrt(r.^2 + r_dash.^2 - 2*(r.*r_dash.*cos(phi_dash - phi)))))./(2.*pi.*sqrt(r.^2 + r_dash.^2 - 2.*(r.*r_dash.*cos(phi_dash - phi)))));
integrated_second_integral = integralN(second_integral,0,a_1,0,2*pi,0,a_2,0,2*pi,'AbsTol',1e-8,'RelTol',1e-3)
Although theoretically the integral of integral3 approach should work, it's somewhat more efficient to do integral2 of integral2. Setting that up is a little more work because integral2 doesn't have an 'ArrayValued' option, so I used the integralN routine I put on the file exchange which does exactly that.
_
>> integrated_second_integral = integralN(second_integral,0,a_1,0,2*pi,0,a_2,0,2*pi,'AbsTol',1e-8,'RelTol',1e-3)
integrated_second_integral =
8.479042788173304e-04 - 2.199114856608778e-08i_
  1 Comment
JJ
JJ on 28 Sep 2014
Mike, thank you for your answer - I have copied the code and module and it works as expected. I now just need to dig into integralN to understand what it's doing :p Thanks again!

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!