Problems in dynamic symsum variables in solving a system of equation

Dear all,

I enounter some problems in using dynamic symsum variables in solving a system of linear equation. I don't know what has gone wrong and don't know how to fix it. My code is as follows:

clear;


syms n
a = sym('a',[1 20])

a = 

b = sym('a_3',[1 20])

b = 

theta_0 = pi/3;
theta_1 = 2*pi/3;
beta = 0.5;
delta = 0.5;
U= delta*exp(i*theta_0) - ((1-delta^2)*exp(i*theta_0/2 + theta_1)/(1-delta*exp(i*theta_1)));

          
eqn1 = beta*a(1) + U*symsum((besselj(1+n ,2)*(-a(n) + beta*b(n))),n,-10,10) == U*b(1)

Error using indexing
Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and function body must be sym expression.

Error in indexing (line 936)
R_tilde = builtin('subsref',L_tilde,Idx);

eqn2 = beta*a(2) + U*symsum((besselj(2+n ,2)*(-a(n) + beta*b(n))),n,-10,10) == U*b(2)
eqn2 = beta*a(3) + U*symsum((besselj(3+n ,2)*(-a(n) + beta*b(n))),n,-10,10) == U*b(2)

Answer 1

[Paul]

Hi Tsz Tsun,

Why not just use regular sum instead of symsum when the goal is to just take the sum of a finite number of terms?

Define a and b, they each have 20 elements

%syms n integer
a = sym('a',[1 20])

a = 

b = sym('a_3',[1 20])

b = 

Proably best to use symbolic constant pi, and symbolic consntants for beta and delta

theta_0 = sym(pi)/3;
theta_1 = 2*sym(pi)/3;
beta = sym(0.5);
delta = sym(0.5);

Used 1i for these expressions, sometimes i is already in use as a variable

U= delta*exp(1i*theta_0) - ((1-delta^2)*exp(1i*theta_0/2 + theta_1)/(1-delta*exp(1i*theta_1)))

U = 

U could probalby be futher simplified.

Here, the original symsum expression was over 21 terms, but a and b only have 20 terms. I assumed they were mathematically indexed from -10 to 9. Also, Matlab uses 1-based indexing, so the indices into a and b are shifted by 11. Use elementwise multiplication, .*, and take the sum. Convert n to sym on the argument to besselj, unless it's desired to express those terms as actual numerical values.

n = -10:9;
eqn1 = beta*a(1) + U*sum(besselj(1+sym(n) ,2).*(-a(n+11) + beta*b(n+11))) == U*b(1)

eqn1 = 

Is that close to the desired result?

%eqn1 = beta*a(1) + U*symsum((besselj(1+n ,2).*(-a(n+11) + beta*b(n+11))),n,-10,9) == U*b(1)

[Tsz Tsun]

Thanks a lot, it is a perfect answer and it is what I need!

[Tsz Tsun]

Now I can proceed to solve my system of equations. For simplicity I run 5 for them (in my actual case I need to run 20 equations to solve for 20 variables). But I get some error when calling the solution, I am not sure what is wrong:

clear;

a = sym('a_1',[1 5])

a = 

b = sym('a_3',[1 5])

b = 

theta_0 = sym(pi)/3;
theta_1 = 2*sym(pi)/3;
beta = sym(0.5);
delta = sym(0.5);
U= delta*exp(1i*theta_0) - ((1-delta^2)*exp(1i*theta_0/2 + theta_1)/(1-delta*exp(1i*theta_1)))

U = 

n = -2:2;
eqn0 = beta*a(1) + U*sum(besselj(0+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(1)

eqn0 = 

eqn1 = beta*a(2) + U*sum(besselj(1+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(2)

eqn1 = 

eqn2 = beta*a(3) + U*sum(besselj(2+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(3)

eqn2 = 

eqn3 = beta*a(4) + U*sum(besselj(3+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(4)

eqn3 = 

eqn4 = beta*a(5) + U*sum(besselj(4+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(5)

eqn4 = 

a(1)= 0;
a(2)= 0; 
a(4)= 0;
a(5)= 0;

S=vpasolve([eqn0, eqn1, eqn2, eqn3,eqn4],...
    [b(1), b(2), b(3), b(4) ,b(5)]);
Sol0 = vpa(S.(b(1)))

Argument to dynamic structure reference must evaluate to a valid field name.

Sol1 = vpa(S.(b(2)))
Sol2 = vpa(S.(b(3)))
Sol3 = vpa(S.(b(4)))
Sol4 = vpa(S.(b(5)))

[Paul]

a = sym('a_1',[1 5]); 
b = sym('a_3',[1 5]);
 
theta_0 = sym(pi)/3;
theta_1 = 2*sym(pi)/3;
beta = sym(0.5);
delta = sym(0.5);
U= delta*exp(1i*theta_0) - ((1-delta^2)*exp(1i*theta_0/2 + theta_1)/(1-delta*exp(1i*theta_1)));

n = -2:2;
eqn0 = beta*a(1) + U*sum(besselj(0+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(1);

eqn1 = beta*a(2) + U*sum(besselj(1+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(2);

eqn2 = beta*a(3) + U*sum(besselj(2+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(3);

eqn3 = beta*a(4) + U*sum(besselj(3+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(4);

eqn4 = beta*a(5) + U*sum(besselj(4+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(5);

The following lines of code didn't do anything insofar as the equations are concerned. If the goal is to set those elements of a to zero in the equations, then need to use subs.

%a(1)= 0;
%a(2)= 0; 
%a(4)= 0;
%a(5)= 0;

subs(eqn0,a([1 2 4 5]),zeros(1,4))

ans = 

No need here write out the elements of b explicitly

%S = vpasolve([eqn0, eqn1, eqn2, eqn3,eqn4],...
%    [b(1), b(2), b(3), b(4) ,b(5)]);

S = vpasolve([eqn0, eqn1, eqn2, eqn3,eqn4],b);

Here, we see that the fields of s are the names of the symbolic variables to which the elements of the matlab variables refer. This solution doesn't reflect the subsitutiond of the zero values for the elements of a.

S

S = struct with fields:
    a_31: - a_11*(1.2329056123032454847802548795858 - 0.13980691173877614830928478085601i) + a_12*(0.76660666261541552495110802368146 - 0.033151883435988784455790965613723i) - a_13*(0.2468087274…
    a_32: a_11*(0.76660666261541552495110802368146 - 0.033151883435988784455790965613723i) - a_12*(0.98915351062819414595173667852215 - 0.12926585899806638779884822874695i) - a_13*(0.97336111976…
    a_33: - a_11*(0.24680872749657166624823350540235 - 0.010673236437883432416009962180343i) - a_12*(0.97336111976951235149822106327399 - 0.042092974086126220156425408379048i) - a_13*(0.93964793…
    a_34: - a_11*(0.82886677465464554064105551840899 - 0.03584432021966388909724197826747i) - a_12*(0.3756823098985887283739206852003 - 0.016246370862770770141491420169485i) - a_13*(0.4404542878…
    a_35: - a_11*(0.52834550533799121559656349320036 - 0.022848286430404755198534440561674i) - a_12*(0.075385836042065797062656340596366 - 0.0032600583468251113173917762224222i) - a_13*(0.157908…

S.a_31

ans = 

%Sol0 = vpa(S.(b(1)))

If we want those substitutions of the zero values for a, then

S = vpasolve(subs([eqn0, eqn1, eqn2, eqn3,eqn4],a([1:2,4:5]),zeros(1,4)),b);
S

S = struct with fields:
    a_31: a_13*(- 0.24680872749657166624823350540235 + 0.010673236437883432416009962180343i)
    a_32: a_13*(- 0.97336111976951235149822106327399 + 0.042092974086126220156425408379048i)
    a_33: a_13*(- 0.93964793756208503536516539027784 + 0.12712499189143282677978044024818i)
    a_34: a_13*(- 0.44045428789206850137011909387728 + 0.019047433218571762240481244454444i)
    a_35: a_13*(- 0.15790840984576462467161045864993 + 0.0068287447162396522877493404774615i)

S.a_31

ans = 

which appears to the same term in the previous solution. Is that expected?

Finallly, I suspect that given the structure of the equations there's probably an easier way to construct them and then manipulate them.

[Tsz Tsun]

Thank you very much for the solution and teaching me to write neat codes!!! This helps perfectly :)

[Paul]

You're weclcome. If, for programmtic reasons for example, you want to reference the fields of S based on the elements of b, the way to do that would be (just creating an S here for an example)

b = sym('a_3',[1 5]);
% create a struct with a fieldname as returned from vpasolve and symbolic
% value in the field
S.a_31 = b(1)*sym(pi)

S = struct with fields:
    a_31: pi*a_31

% retreive
S.(string(b(1)))

ans = 

[Tsz Tsun]

Thanks a lot, it really helps !

[Paul]

a = sym('a_1',[1 5]); 
b = sym('a_3',[1 5]);
 
theta_0 = sym(pi)/3;
theta_1 = 2*sym(pi)/3;
beta = sym(0.5);
delta = sym(0.5);
U= delta*exp(1i*theta_0) - ((1-delta^2)*exp(1i*theta_0/2 + theta_1)/(1-delta*exp(1i*theta_1)));
n = -2:2;

Original formualation

eqn0 = beta*a(1) + U*sum(besselj(0+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(1);
eqn1 = beta*a(2) + U*sum(besselj(1+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(2);
eqn2 = beta*a(3) + U*sum(besselj(2+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(3);
eqn3 = beta*a(4) + U*sum(besselj(3+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(4);
eqn4 = beta*a(5) + U*sum(besselj(4+sym(n) ,2).*(-a(n+3) + beta*b(n+3))) == U*b(5);
S = vpasolve([eqn0, eqn1, eqn2, eqn3,eqn4],b);

Vectorized formulation, one line

eqns = beta*a(:) + U*sum(besselj((0:4).'+sym(n) ,2).*(-a(n+3) + beta*b(n+3)),2) == U*b(:);

Solve and compare

Snew = vpasolve(eqns,b);
isequal(Snew,S)

ans = logical
   1

The function equationsToMatrix also may be of interest.

[Tsz Tsun]

The Vectorized formulation seems to be quite advanced to me, it is very elegant, thanks so much for introducing this techinque to me! But I don't quite get the .'+sym(n) syntax, what is the need for ' before + ?

[Paul]

The operator .' is shorthand for transpose, .' (which is, in general, not the same as ctranspose, '). The code is taking advantage of "binary singleton expansion," see here: Compatible Array Sizes.

(0:4).' + sym(n)

is adding a column vector and row vector to yield a matrix, which in turn is used to compute besselj at each element of thta matrix.

That matrix is then elementwise muliplied with .* with (-a(n+3) + beta*b(n+3)), which is a row vector. That row vector is "expanded" by copying that row as many time as needed to form a matrix that has as many rows as the matrix returned from besselj. That result is then element-wise multiplied with the besselj matrix. The resulting matrix is summed across the columns to yield a column vector.

If still not clear, break up that line into multiple steps, one for each operation.

[Tsz Tsun]

Thanks for the explanation ! It is very clear!

Answer 2

[Walter Roberson]

A symbolic variable or expression or symfun or symmatrix can never be used as an index.

I cannot show you the corrected code because you are also trying to use n as negative indices -- you want to use a(n) with n from -10 to 10, but there is no possibility of a(-10)

You could construct

N = -10:10;

A = sym("a_" + regexprep(string(N), '-', '_'))
B = sym("b_" + regexprep(string(N), '-', '_'))

an = @(N) A(N+(length(A)+1)/2)
bn = @(N) B(N+(length(B)+1)/2)

and then you would be able to use things like

sum(besselj(2+N).*(-a(N) + beta*b(N)))

[Tsz Tsun]

Thanks a lot!

[Walter Roberson]

I corrected the an and bn definitions -- had a negative that should not have been there.

[Tsz Tsun]

should it be sum(besselj(1+N).*(-a(N) + beta*b(N))) also instead?

[Walter Roberson]

No, before you were trying to use besselj(2+n,2) with n starting from -10, so when N is -10 to +10 then you would use besselj(2+N,2) not besselj(1+N,2)

I did, however, miss the ,2 parameter for besselj

[Tsz Tsun]

I see, thanks for the clarification !