Unable to use quadgk for integration

Question

RB 2017-5-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/342499-unable-to-use-quadgk-for-integration

评论： RB 2017-5-30

I have been struggling to use quadgk for computing the integral ogf a function.

The code and the error are given below:

if true
  n=35; 
T=15; 
g=.111;
K=1/g;
d=2; 
rbar=0.04;
mubar=0; 
R=[1, 0; 0, 0];
M=[0, 0; 0, 1];
X012=0.002;
X0=[0.01, X012; X012, 0.001];
H=[-0.5, 0.4; 0.007, -0.008];
Q=[0.06 -0.0006; -0.06, 0.006];
Scheck7=0;
beta=3;
SIGMA=[0.006811791307233, -0.000407806990090; -0.000407806990090, 0.00039291243623];
PSICURL=[0.011526035149236, 0.758273970303934; 0.013935191202735, 0.955423605940771];   
S0i=zeros(n,1);
S0i2=zeros(n,1);
Yn0=0; %Used in the first lower bound
Yi=zeros(n,1);
SPOW=0;
SPHI=0;
SPSI=0;
% Parameters used for the new approach
% To define the Matrices PSIi in an array
delta=0.75; %The damping factor
PSIi = cell(1, n);
PHIi=zeros(n,1);
% Specification of new MATRICES Aij's
ANEW=expm([T.*H,T.*(2*(Q'*Q));T.*(R+M),T.*(-(H'))]);
A11=ANEW(1:2,1:2); %That is good
A21=ANEW(3:4,1:2);
A12=ANEW(1:2,3:4);%intersection
A22=ANEW(3:4,3:4);
% Computationof psi(T) and phi(T) (Otherwise with varying i;run in a loop)
C22=inv(A22);
PSIT=C22*A21;
PHIT=beta*(log(det(A22))+T*trace(H'))/2;
% Computation of SZCB's Pcurl(0,T)
C=trace(PSIT*X0);
SZCB=exp(-(rbar+mubar)*T)*exp(-PHIT-C)
SIGMAi=inv(SIGMA);
THETA1=SIGMAi*(PSICURL'*X0*PSICURL);
for i=2:n;
       APHNEW=expm([(i-1).*H,(i-1).*(2*(Q'*Q));(i-1).*(R+M),(i-1).*(-(H'))]);
       APHNEW11=APHNEW(1:2,1:2); %That is good
       APHNEW21=APHNEW(3:4,1:2);
       APHNEW12=APHNEW(1:2,3:4);
       APHNEW22=APHNEW(3:4,3:4);
       % Computation of PSIi and PHIi  
       BPHNEW22=inv(APHNEW22);
       PSIi{i}=APHNEW22\APHNEW21;
       M7=PSIi{i};
       SPSI=SPSI+PSIi{i};
       PHIi(i)=beta*(log(det(APHNEW22))+(i-1)*trace(H'))/2; %comparing with mcir the trace takes care of the loop
       S0i(i)=exp(-((rbar+mubar)*(i-1)+PHIi(i)));
       S0i2(i)=(-((rbar+mubar)*(i-1)+PHIi(i))); %S0i without the exponent
       Yn0=Yn0+S0i2(i);
end
SIGMAi=inv(SIGMA);
THETA1=SIGMAi*(PSICURL'*X0*PSICURL);
fun0=@(Gamma) exp(trace(1i.*(THETA1*inv(eye(2)-2i.*SIGMA*Gamma)*SIGMA*Gamma)))/((det(eye(2)-2i.*SIGMA*Gamma))^(beta/2));
fun4 = @(Gamma1) exp(-(1i*Gamma1)*log(K-1))*exp((1i*Gamma1+(delta+1))*Yn0)*(fun0((-(Gamma1-1i*(delta+1))).*SPSI))/(delta^2+delta-Gamma1^2+(1i*Gamma1*(2*delta+1)));
qty = quadgk(fun4,0,inf)
qty1= exp(-(delta*log(K-1)))*qty/pi
LB1=g*SZCB*qty1
format 'long'
end

The error I get is: if true

Error using * Inner matrix dimensions must agree. Error in lkag7>@(Gamma1)exp(-(1i*Gamma1)*log(K-1))*exp((1i*Gamma1+(delta+1))*Yn0)*(fun0((-(Gamma1-1i*(delta+1))).*SPSI))/(delta^2+delta-Gamma1^2+(1i*Gamma1*(2*delta+1))) Error in quadgk/evalFun (line 330) fx = FUN(x); Error in quadgk/f2 (line 361) [y,too_close] = evalFun(t2t); Error in quadgk/vadapt (line 249) [fx,too_close] = f(x); Error in quadgk (line 196) [q,errbnd] = vadapt(@f2,interval); Error in weblkag7 (line 67) qty = quadgk(fun4,0,inf) end

I would be grateful if anyone can help.

Thanks