complex integral problem

Hi all,

I have to solve following complex integral with 4 variables (s1x,s1y,s2x,s2y) numerically in Matlab GUI. But i i have error and couldn’t find the error till now. And codes are as follows too.

Could anyone show me the error?

Thanks in advance.

??? Error using ==> quadgk at 140

Parameters must be a character.

Error in ==> TransCalc>ApproveBut_Callback at 784

      TransCoeff{m}=quadgk(y{m},handles.S1xMinVal,handles.S1xMaxVal,handles.S1yMinVal,handles.S1yMaxVal,handles.S2xMinVal,handles.S2xMaxVal,handles.S2yMinVal,handles.S2yMaxVal)

Error in ==> gui_mainfcn at 96

          feval(varargin{:});

Error in ==> TransCalc at 42

      gui_mainfcn(gui_State, varargin{:});

Error in ==> @(hObject,eventdata)TransCalc('ApproveBut_Callback',hObject,eventdata,guidata(hObject))

??? Error while evaluating uicontrol Callback

for d=floor(10000/handles.WaveUp):1:ceil(10000/handles.WaveLow)

a(m)=d;

k(m)=2*pi/((10000*(1e-6)/d));

lamda(m)=10000*(1e-6)/d;

n_alfa(m)=0.5*gamma(handles.AlfaPowLaw)*((2*pi)^(-11/6+(handles.AlfaPowLaw)/2)).*((lamda(m).*(handles.DistLow))^(11/6-(handles.AlfaPowLaw)/2));

d_alfa=gamma(1-handles.AlfaPowLaw/2)*((gamma(handles.AlfaPowLaw/2))^2)*gamma(handles.AlfaPowLaw-1)*cos(pi*handles.AlfaPowLaw/2)*sin(pi*handles.AlfaPowLaw/4);

cn2tilda(m)=(n_alfa(m)/d_alfa)*handles.StrConst; 
set(handles.StrucConstNonKolmNum,'String',num2str(d_alfa));

ro_0(m)=(((-2)^(1-handles.AlfaPowLaw))*gamma((2-handles.AlfaPowLaw)/2)*gamma(handles.AlfaPowLaw-1)*cos(0.5*pi*handles.AlfaPowLaw)*((k(m))^2)*handles.DistLow*((cn2tilda(m))^2)/gamma(handles.AlfaPowLaw/2))^(-1/(handles.AlfaPowLaw-2));

f1=exp((-1/((handles.BeamRad)^2))*(((s1x)^2)+((s1y)^2)+((s2x)^2)+((s2y)^2)));

f2{m}=exp((1i*(k(m)/(2*handles.DistLow)))*((((s1x)^2)-(2*s1x*handles.PxAxisValue))+(((s1y)^2)-(2*s1y*handles.PyAxisValue))+((2*s2x*handles.PxAxisValue)-((s2x)^2))+((2*s2y*handles.PyAxisValue)-((s2y)^2))))

f3{m}=exp((-1/((ro_0(m))^(handles.AlfaPowLaw-2)))*((((s1x-s2x)^2)+((s1y-s2y)^2))^((handles.AlfaPowLaw-2)/2)));

f4(m)=((handles.DistLow*handles.DistLow+k(m)*k(m)*((handles.BeamRad)^4)))/((4*pi*pi*handles.DistLow*handles.DistLow*((handles.BeamRad)^4))*((handles.PxAxisValue*handles.PxAxisValue+handles.PyAxisValue*handles.PyAxisValue)*((-k(m)*k(m)*(handles.BeamRad)*(handles.BeamRad))/(handles.DistLow*handles.DistLow+k(m)*k(m)*((handles.BeamRad)^4)))));

y{m} = @(s1x,s1y,s2x,s2y) (f1.*f2{m}.*f3{m})

TransCoeff{m}=quadgk(y{m},handles.S1xMinVal,handles.S1xMaxVal,handles.S1yMinVal,handles.S1yMaxVal,handles.S2xMinVal,handles.S2xMaxVal,handles.S2yMinVal,handles.S2yMaxVal)

m=m+1;

end

Answer 1

[Walter Roberson]

You cannot pass variables to quadgk that way. See instead http://www.mathworks.com/help/techdoc/math/bsgprpq-5.html

[Walter Roberson]

quadgk is *not* a routine for integrating over 4 variables: the "quad" refers to adaptive quadrature.

There are techniques for doing multiple integrals, starting from the lower-dimensional routines such as triplequad . I do not recall the method at the moment. If I recall correctly, Mike Hosa has shown a general method, but I do not find it right now.

Answer 2

[Yalcin]

Thanks for your answer. I changed 4-d integration as follow. But now i have following error. what should i do?

for d=floor(10000/handles.WaveUp):1:ceil(10000/handles.WaveLow) a(m)=d; k(m)=2*pi/((10000*(1e-6)/d)); lamda(m)=10000*(1e-6)/d; n_alfa(m)=0.5*gamma(handles.AlfaPowLaw)*((2*pi)^(-11/6+(handles.AlfaPowLaw)/2)).*((lamda(m).*(handles.DistLow))^(11/6-(handles.AlfaPowLaw)/2)); d_alfa=gamma(1-handles.AlfaPowLaw/2)*((gamma(handles.AlfaPowLaw/2))^2)*gamma(handles.AlfaPowLaw-1)*cos(pi*handles.AlfaPowLaw/2)*sin(pi*handles.AlfaPowLaw/4); cn2tilda(m)=(n_alfa(m)/d_alfa)*handles.StrConst;

    ro_0(m)=(((-2)^(1-handles.AlfaPowLaw))*gamma((2-handles.AlfaPowLaw)/2)*gamma(handles.AlfaPowLaw-1)*cos(0.5*pi*handles.AlfaPowLaw)*((k(m))^2)*handles.DistLow*((cn2tilda(m))^2)/gamma(handles.AlfaPowLaw/2))^(-1/(handles.AlfaPowLaw-2));

    f1=@(s1x,s1y,s2x,s2y)exp((-1/((handles.BeamRad)^2))*(((s1x)^2)+((s1y)^2)+((s2x)^2)+((s2y)^2)));

    f2(m)=@(s1x,s1y,s2x,s2y)exp((1i*(k(m)/(2*handles.DistLow)))*((((s1x)^2)-(2*s1x*handles.PxAxisValue))+(((s1y)^2)-(2*s1y*handles.PyAxisValue))+((2*s2x*handles.PxAxisValue)-((s2x)^2))+((2*s2y*handles.PyAxisValue)-((s2y)^2))));

    f3(m)=@(s1x,s1y,s2x,s2y)exp((-1/((ro_0(m))^(handles.AlfaPowLaw-2)))*((((s1x-s2x)^2)+((s1y-s2y)^2))^((handles.AlfaPowLaw-2)/2)));

    f4(m)=((handles.DistLow*handles.DistLow+k(m)*k(m)*((handles.BeamRad)^4)))/((4*pi*pi*handles.DistLow*handles.DistLow*((handles.BeamRad)^4))*((handles.PxAxisValue*handles.PxAxisValue+handles.PyAxisValue*handles.PyAxisValue)*((-k(m)*k(m)*(handles.BeamRad)*(handles.BeamRad))/(handles.DistLow*handles.DistLow+k(m)*k(m)*((handles.BeamRad)^4)))));

    y(m) = @(s1x,s1y,s2x,s2y)(f4(m)*f1*f2(m)*f3(m));

    tol=1e4;

    %f5(m) = @(s1x,s1y,s2x) quadv(@(s2y) y(m),handles.S2yMinVal,handles.S2yMaxVal,tol);
    %TransCoeff = triplequad(f5(m), handles.S1xMinVal,handles.S1xMaxVal,handles.S1yMinVal,handles.S1yMaxVal,handles.S2xMinVal,handles.S2xMaxVal, tol);
   TransCoeff = triplequad((@(s1x,s1y,s2x) quadv(@(s2y) y(m),handles.S2yMinVal,handles.S2yMaxVal,tol)), handles.S1xMinVal,handles.S1xMaxVal,handles.S1yMinVal,handles.S1yMaxVal,handles.S2xMinVal,handles.S2xMaxVal,tol);

    m=m+1;
end

??? Undefined function or method 'mtimes' for input arguments of type 'function_handle'.

Error in ==> TransCalc>@(s1x,s1y,s2x,s2y)(f4(m)*f1*f2(m)*f3(m)) at 804 y(m) = @(s1x,s1y,s2x,s2y)(f4(m)*f1*f2(m)*f3(m));

Error in ==> TransCalc>@(s2y)y(m) at 811 TransCoeff = triplequad((@(s1x,s1y,s2x) quadv(@(s2y) y(m),handles.S2yMinVal,handles.S2yMaxVal,tol)), handles.S1xMinVal,handles.S1xMaxVal,handles.S1yMinVal,handles.S1yMaxVal,handles.S2xMinVal,hand Error in ==> quadv at 57 y{j} = feval(f, x(j), varargin{:}); %#ok<AGROW>

Error in ==> TransCalc>@(s1x,s1y,s2x)quadv(@(s2y)y(m),handles.S2yMinVal,handles.S2yMaxVal,tol) at 811 TransCoeff = triplequad((@(s1x,s1y,s2x) quadv(@(s2y) y(m),handles.S2yMinVal,handles.S2yMaxVal,tol)), handles.S1xMinVal,handles.S1xMaxVal,handles.S1yMinVal,handles.S1yMaxVal,handles.S2xMinVal,hand Error in ==> triplequad>innerintegral at 59 fcl = intfcn(xmin, y(1), z(1), varargin{:});

Error in ==> dblquad>innerintegral at 74 fcl = intfcn(xmin, y(1), varargin{:}); %evaluate only to get the class below

Error in ==> quad at 77 y = f(x, varargin{:});

Error in ==> dblquad at 54 Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

Error in ==> triplequad at 47 Q = dblquad(@innerintegral, ymin, ymax, zmin, zmax, tol, quadf, intfcn, ...

Error in ==> TransCalc>ApproveBut_Callback at 811 TransCoeff = triplequad((@(s1x,s1y,s2x) quadv(@(s2y) y(m),handles.S2yMinVal,handles.S2yMaxVal,tol)), handles.S1xMinVal,handles.S1xMaxVal,handles.S1yMinVal,handles.S1yMaxVal,handles.S2xMinVal,hand Error in ==> gui_mainfcn at 96 feval(varargin{:});

Error in ==> TransCalc at 42 gui_mainfcn(gui_State, varargin{:});

Error in ==> @(hObject,eventdata)TransCalc('ApproveBut_Callback',hObject,eventdata,guidata(hObject))

??? Error while evaluating uicontrol Callback