Error using sym/subsindex (line 737)

Hello!

I am trying to execute the following code:

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; 
Yi=zeros(n,1);
S1MC=zeros(n,1);
POW=zeros(n,1);
SPOW=0;
SPHI=0;
SPSI=0;
delta=0.75; 
PSIi = cell(1, n);
PHIi=zeros(n,1);
% Prameters for the symbolic sum
syms Gamma1 ik
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; 
         S0i(i)=exp(-((rbar+mubar)*(i-1)+PHIi(i)));
         S0i2(i)=(-((rbar+mubar)*(i-1)+PHIi(i))); 
         Yn0=Yn0+S0i2(i);
end
% i stands for iota
fun0=@(Gamma) exp(trace(1i*(THETA1*inv(eye(2)-2i*SIGMA*Gamma)*SIGMA*Gamma)))/((det(eye(2)-2i*SIGMA*Gamma))^(beta/2));
fun3 = matlabFunction((exp((1i*Gamma1+delta)*Yn0)*(symsum(S0i(ik)*fun0(1i*PSIi{ik}-(Gamma1-1i*delta)*SPSI), ik, 2, n)-(K-1)*fun0(-(Gamma1-1i*delta)*SPSI)))/(1i+Gamma1));
a=fun3(1)
q = quadgk(fun3,0,Inf)
q1= exp(-(delta*x1))*max(q/pi, 0)
x = fminbnd(q1,-Inf,+Inf)

*Please see attached file***

and get the following error:

Error using sym/subsindex (line 737)
Invalid indexing or function definition. When defining a function, ensure that the arguments are symbolic variables and
the body of the function is a SYM expression. When indexing, the input must be numeric, logical, or ':'.
Error in fweb7 (line 51)
fun3 = matlabFunction((exp((1i*Gamma1+delta)*Yn0)*(symsum(S0i(ik)*fun0(1i*PSIi{ik}-(Gamma1-1i*delta)*SPSI), ik, 2,
n)-(K-1)*fun0(-(Gamma1-1i*delta)*SPSI)))/(1i+Gamma1));

Can anyone please help?

Thanks in advance.

Answer 1

[Walter Roberson]

Your fun3 line includes

symsum(S0i(ik)*fun0(1i*PSIi{ik}-(Gamma1-1i*delta)*SPSI), ik, 2, n)

Inside the symsum, you try to use the symbol ik to index PSIi .

In MATLAB, it is never possible to use a symbol to index something -- not possible for vectors, not possible for cell arrays.

I have attached a re-worked version. It is not possible to do the integration using quadqk because at the point of the integration you still have an unknown parameter, THETA1, that you will be one of the two parameters you are minimizing over.

I would note, though, that your exp(-(delta*x1)) and your max(q,0) are separable, so you can reduce the work required for minimization. exp(-(delta*x1)) is strictly >= 0 for real-valued x1. max(q,0) is >= 0 for all q. The product must then be >= 0. You could achieve the minimum of the product if either value were 0.

q looks complicated to compute and solve for a root, but the root of exp(-(delta*x1)) is easy to compute: for positive delta, it occurs when x1 is infinity, giving exp(-infinity) which is 0. Your fminbnd includes the entire range -inf to +inf so +inf is a valid value for x1.

Therefore, by examination we can solve the entire problem as x = +infinity, without ever having computed q .

The only exception would be if we could prove that q is everywhere +infinity, over all THETA1, as in that case you would be working with zero times infinity which is undefined. Therefore if we could prove that q is +infinity for all THETA1, we would have to prohibit x1 = infinity, and so exp(-(delta*x1)) would be positive; in that case, since the positive value is being multiplied by infinity (by hypothesis), it does not matter what value of x1 we use (other that +infinity) as all of them would minimize the product to the same value, +infinity.

You should probably proceed from here to prove that at least one value of q is finite for some THETA1, to make the overall solution x1 = +infinity valid.

[RB]

Thanks for the answer. I still have a few problems:

1)THETA1 is known. I missed on two lines while I was creating the web code.

if***
 SIGMAi=inv(SIGMA);
THETA1=SIGMAi*(PSICURL'*X0*PSICURL);  
end

2)In fact x1 is unknown. x1 is a scalar. I gave an arbitrary value of 1 to x1 in the previous code to see if quadgk works.

P.S. Attached question file.

3)I use your code with modifications by replacing THETA1 by x1 but do not get any answer.

P.S. Attached Code.

Thanks.

[Walter Roberson]

Adjusted version attached

You will find that the answer it gets is the same as the random guess it takes for the minimization. This is because q comes out complex, so max(q/pi,0) is complex (and the same as q/pi no matter what the signs of the components of q are), so you are attempting to do a minimization of a complex function, and that returns immediately.

If q had been real-valued instead of complex then the result would be as I described before, that you could skip the entire calculation of q and instead use x1 = infinity .

[RB]

Hello! Thanks once again.

Is there any other option as 'vpaintegral' is not being recognized.

I have MATLAB2015b.

I get the following error.

if true
  % code
end
Undefined function or variable 'vpaintegral'.
Error in THEFINAL7 (line 59)
q_formula = vpaintegral(fun3_formula, Gamma1, 0, inf, 'RelTol', 1E-30);
Thanks in advance.

[Walter Roberson]

fun3 = matlabFunction(fun3_formula);
q = integral(fun3, 0, inf);

[RB]

Thanks once again. It works giving me a complex answer. But I still have one doubt. When you define your fun3_formula

if true
  % code
fun3_formula = ((exp(-(1i*Gamma1)*x1))*exp((1i*Gamma1+delta)*Yn0)*(sum3a-(K-1)*fun0(-(Gamma1-1i*delta)*SPSI)))/(1i+Gamma1);

you do not consider

if true
  % code
exp(-(1i*Gamma1)*x1.

Can you elaborate the reason for this? because if I include this with

if true
  % code
syms x1
end

I get the following error.

if true
  Error using
Not enough input arguments.
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 84)
      [q,errbnd] = vadapt(@AToInfInvTransform,interval);
Error in integral (line 89)
Q = integralCalc(fun,a,b,opstruct);
Error in THEFINAL7 (line 81)
q = integral(fun3, 0, inf);  end

I am expecting a real answer. Thanks in advance.

[Walter Roberson]

Your original code had

fun3 = matlabFunction((exp((1i*Gamma1+delta)*Yn0)*(symsum(S0i(ik)*fun0(1i*PSIi{ik}-(Gamma1-1i*delta)*SPSI), ik, 2, n)-(K-1)*fun0(-(Gamma1-1i*delta)*SPSI)))/(1i+Gamma1));
a=fun3(1)
q = quadgk(fun3,0,Inf)
q1= exp(-(delta*x1))*max(q/pi, 0)

so your fun3 does not include x1. You do not integrate with respect to x1: you calculate q as a numeric constant and minimize what is effectively

if q <= 0
  q1 = @(x1) zeros(size(x1));
else
  q1 = @(x1) exp(-(delta*x1))*q/pi;
end

but since q is constant minimizing q1 is the same as

q/pi * min( exp(-(delta*x1)) over x1 = -inf to inf)

but as I discussed before, that will come at x1 = +inf, giving you q/pi * 0 for an overall minimum of 0, with no point in doing any of the previous calculations.

[RB]

Thanks for all your help. I completely accept your answer.

Does this mean that with x1 included it is not possible to obtain a maximum?

P.S. question7.pdf file attached above.

Thanks again.