How to use loop with vpasolve ??

In the following code, I want to use loop with vpasolve. For different values of alpha, I want to solve the equation to get different values of r.

syms alpha mio B r
mio=0.6;
B=2;
alpha=1:0.5:6;
alpha_length=length(alpha);
          for i=1:alpha_length
                    bast(i)=(r(i)/0.5).^(alpha(i)*(mio-1));
                    mqam_part1(i)=3*B*((sqrt(3)/2).^(alpha(i)*mio));
                    mqam_part2(i)=((0.5*sqrt(3)).^(-alpha(i)))+((1.5*sqrt(3)).^(-alpha(i)));
                    mqam_part3(i)=3*(((r(i)/0.5)).^(alpha(i)*mio));
                    mqam_part4(i)=((sqrt(3)-(r(i)/0.5)).^-alpha(i));
                    mqam_part5(i)=((2*sqrt(3))-(r(i)/0.5)).^-alpha(i);
                    mqam_part6(i)=mqam_part4(i)+mqam_part5(i);
                    mqam_part7(i)=2*((3-(r(i)/0.5)).^-alpha(i));
                    mqam_part8(i)=6*(((r(i)/0.5)).^(alpha(i)*mio));
                    mqam_part9(i)=6*B*(2.^-alpha(i));
                    eqn_LHS(i)=bast(i)/(mqam_part1(i)+mqam_part2(i))+(mqam_part3(i)*(mqam_part6(i)+mqam_part7(i)));
                    eqn_RHS(i)=B/((mqam_part8(i)*mqam_part6(i))+mqam_part9(i));
                    eqn1(i)=eqn_LHS(i)==eqn_RHS(i);
                    sol_positive(i) = vpasolve(eqn1(i),r(i),[0 Inf]);
          end

But, After running this code, it shows the following error.

Error using subsref
Index exceeds matrix dimensions.

Error in sym/subsref (line 771)
            R_tilde = builtin('subsref',L_tilde,Idx);

Error in tetra_v3 (line 10)
                    bast(i)=(r(i)/0.5).^(alpha(i)*(mio-1));

So, my question is: What is the cause of these errors and how to solve them ?? and How to use loop with vpasolve ??

[Image Analyst]

Original question by Erman:

In the following code, I want to use loop with vpasolve. For different values of alpha, I want to solve the equation to get different values of r.

syms alpha mio B r
mio=0.6;
B=2;
alpha=1:0.5:6;
alpha_length=length(alpha);
          for i=1:alpha_length
                    bast(i)=(r(i)/0.5).^(alpha(i)*(mio-1));
                    mqam_part1(i)=3*B*((sqrt(3)/2).^(alpha(i)*mio));
                    mqam_part2(i)=((0.5*sqrt(3)).^(-alpha(i)))+((1.5*sqrt(3)).^(-alpha(i)));
                    mqam_part3(i)=3*(((r(i)/0.5)).^(alpha(i)*mio));
                    mqam_part4(i)=((sqrt(3)-(r(i)/0.5)).^-alpha(i));
                    mqam_part5(i)=((2*sqrt(3))-(r(i)/0.5)).^-alpha(i);
                    mqam_part6(i)=mqam_part4(i)+mqam_part5(i);
                    mqam_part7(i)=2*((3-(r(i)/0.5)).^-alpha(i));
                    mqam_part8(i)=6*(((r(i)/0.5)).^(alpha(i)*mio));
                    mqam_part9(i)=6*B*(2.^-alpha(i));
                    eqn_LHS(i)=bast(i)/(mqam_part1(i)+mqam_part2(i))+(mqam_part3(i)*(mqam_part6(i)+mqam_part7(i)));
                    eqn_RHS(i)=B/((mqam_part8(i)*mqam_part6(i))+mqam_part9(i));
                    eqn1(i)=eqn_LHS(i)==eqn_RHS(i);
                    sol_positive(i) = vpasolve(eqn1(i),r(i),[0 Inf]);
          end

But, After running this code, it shows the following error.

Error using subsref
Index exceeds matrix dimensions.

Error in sym/subsref (line 771)
            R_tilde = builtin('subsref',L_tilde,Idx);

Error in tetra_v3 (line 10)
                    bast(i)=(r(i)/0.5).^(alpha(i)*(mio-1));

So, my question is: What is the cause of these errors and how to solve them ?? and How to use loop with vpasolve ??

Answer 1

[John D'Errico]

What is the cause? READ THE ERROR MESSAGE.

"Index exceeds matrix dimensions."

What are you indexing?

bast(i)=(r(i)/0.5).^(alpha(i)*(mio-1));

What indexing is involved here? alpha seems to be a vector.

How about r?

r is a scalar. When i is greater than 1, r(i) will return an error, because r IS A SCALAR. What is the second element of a scalar?

[Eman S]

Thanks for reply. How can I solve this error and how can I solve this equation using loop ??

[John D'Errico]

How can you solve it? r is not a vector. You can't solve it, at least not with the code you have written. r has only one value.

If r has only one value, why are you trying to use multiple values for r? Why do you need to index r at all?

[Eman S]

I want to solve an equation for multiple times to get the value of r by varying the value of alpha in each time. The values of alpha are 1:0.5:6. With each value of alpha, i want to solve the equation to get the value of r.

[John D'Errico]

You are saving the solution in sol_positive(i), NOT in r(i). There is no need to index r. r is just the symbolic variable you are solving for in the equation.

[Eman S]

OK. I modified the code as the following. I removed indexing from r.

syms  mio B r
mio=0.6;
B=2;
alpha=1:0.5:6;
alpha_length=length(alpha);
          for i=1:alpha_length
                    bast(i)=(r/0.5).^(alpha(i)*(mio-1));
                    mqam_part1(i)=3*B*((sqrt(3)/2).^(alpha(i)*mio));
                    mqam_part2(i)=((0.5*sqrt(3)).^(-alpha(i)))+((1.5*sqrt(3)).^(-alpha(i)));
                    mqam_part3(i)=3*(((r/0.5)).^(alpha(i)*mio));
                    mqam_part4(i)=((sqrt(3)-(r/0.5)).^-alpha(i));
                    mqam_part5(i)=((2*sqrt(3))-(r/0.5)).^-alpha(i);
                    mqam_part6(i)=mqam_part4(i)+mqam_part5(i);
                    mqam_part7(i)=2*((3-(r/0.5)).^-alpha(i));
                    mqam_part8(i)=6*(((r/0.5)).^(alpha(i)*mio));
                    mqam_part9(i)=6*B*(2.^-alpha(i));
                    eqn_LHS(i)=bast(i)/(mqam_part1(i)+mqam_part2(i))+(mqam_part3(i)*(mqam_part6(i)+mqam_part7(i)));
                    eqn_RHS(i)=B/((mqam_part8(i)*mqam_part6(i))+mqam_part9(i));
                    eqn1(i)=eqn_LHS(i)==eqn_RHS(i);
                    sol_positive(i) = vpasolve(eqn1(i),r,[0 Inf]);
          end

After running this code, it gives the following errors:

Error using subsasgn
In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in sym/privsubsasgn (line 997)
                L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

Error in sym/subsasgn (line 834)
            C = privsubsasgn(L,R,inds{:});
Error in tetra_v3 (line 23)
                    sol_positive(i) = vpasolve(eqn1(i),r,[0 Inf]);

Answer 2

[Walter Roberson]

As we explored earlier, your system works out to be a polynomial and vpasolve() is going to return a list of all the solutions under the constraint you give, [0 inf]. You are trying to store that vector into a single location sol_positive(i) .

If you were certain that there would always be the same number of results, you could store to sol_positive(i,:) instead, but I think you would be better off assuming that the number of positive roots might change, so I would recommend assigning to sol_positive{i} . Indeed, my test shows that most of your equations have no solution in that range.

[Walter Roberson]

Correction: this is a different system that is not polynomial.

However, your system does not happen to have solutions at the alpha that end in 0.5 . There are solutions with fractional alpha, but those solutions do not happen be real valued for alpha ending in 1/2 . For example for r = 1/5 there is a solution of alpha about 2.4

Also, because it is not polynomial, vpasolve() is only finding one solution. For alpha = 1, there are three positive solutions in the range 1 to 1.7