Unable to use quadgk for integration

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

Answer 1

[Torsten]

In fun0 and fun4, replace all "*" , "/" and "^" by ".*", "./" and ".^" .

Best wishes

Torsten.

[RB]

Thanks for the answer. But since matrices are involved and .* would mean element wise operations. Would it be correct?

Thanks in advance.

Even if I change to .* ./ and .^ it gives the same error.

I use

if true
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)));
end

The same problem appears.

[Torsten]

"quadgk" will supply a vector for 'gamma' and expects a vector of function values to be returned by 'fun4'.

Just test whether your function works by defining a vector for 'gamma' and setting y = fun4(gamma).

Best wishes

Torsten.

[RB]

Thanks. Gamma is a vector and gamma1 is a scalar. I test with gamma1 = 2 and get the following

if true
  a=fun4(2)
end

if true
  a =

      1.749250335010293e-33 - 5.565461789832223e-34i
  end

[Torsten]

I repeat:

"fun4" will be called with a vector x, and you will have to supply the corresponding vector of function values y = fun4(x).

Best wishes

Torsten.

[RB]

Thanks again.

However I do not have vector x. gamma and gamma1 are unknown. I want to integrate w.r.t. gamma1.

If I merge fun4 and fun0 as below, I still get the error

if true
  fun4 = @(Gamma1) exp(-(1i.*Gamma1).*log(K-1)).*exp((1i.*Gamma1+(delta+1)).*Yn0).*(exp(trace(1i.*(THETA1.*inv(eye(2)-2i.*SIGMA.*((-(Gamma1-1i.*(delta+1))).*SPSI)).*SIGMA.*((-(Gamma1-1i.*(delta+1))).*SPSI))))./((det(eye(2)-2i.*SIGMA.*((-(Gamma1-1i.*(delta+1))).*SPSI))).^(beta./2)))./(delta.^2+delta-Gamma1.^2+(1i*Gamma1.*(2.*delta+1)));  end

[Torsten]

As written in your code, you try to integrate "fun4" over [0:Inf).

For this purpose, "quadgk" calls "fun4" with a vector 'gamma1' taken from the interval [0:Inf) and expects you to supply fun4(gamma1) for this vector 'gamma1'.

I don't know where your functions "fun0" and "fun4" stem from, but that's the way "quadgk" will work.

Best wishes

Torsten.

[RB]

Thanks so much for replying. I am testing a lot but I have a querry that whenever I try to include the trace in exponential, I start getting the same error. If you see, when I merge the two functions, I will need a scalar gamma1.

[Torsten]

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.

[RB]

Thanks again. Can you have a look at the attached code?