Unable to use quadgk for integration
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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
0 comentarios
Respuestas (1)
Torsten
el 30 de Mayo de 2017
In fun0 and fun4, replace all "*" , "/" and "^" by ".*", "./" and ".^" .
Best wishes
Torsten.
9 comentarios
Torsten
el 30 de Mayo de 2017
Then write a function routine for "fun4" in which you loop over the elements of gamma1. This way, you only work with scalars.
Best wishes
Torsten.
Ver también
Categorías
Más información sobre Gamma Functions en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!