How to use a symbolic variable to refer to a specific cell on an array

I am trying to plot this function for ρ=constant(380) and z=constant(12) (so as a function of φ) but when I run the code I have a problem because i use a double symsum to refer to the Smn values of an array. I know that a symbolic variable cannot be used as an index of an array but I don't know how to overcome this problem. Jm is the mth Bessel waveform, Smn is the nth root of the mth Bessel waveform.The file that generates the J array is the BesselJ.m and the file that generates the roots of the mth bessel function is the Smn.n along with his function bz.m.

The error code is: Error using sym/subsindex (line 855)

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 Untitled11 (line 10)

b=symsum(symsum(vpaintegral(J(2*m-1,:)*380/Ro,r,[0

Ro]).*interp1(v,J(2*m-1,:),s(2*m-1,2*n-1)*380/Ro).*cos((2*n-1)*pi*phi/tp)*(exp(s(2*m-1,2*n-1)*z/Ro)+exp((2*l-z)*s(2*m-1,2*n-1)/Ro))./((2*n-1)*pi*interp1(v,J(2*m,:).^2,s(2*m-1,2*n-1))),n,1,10),m,1,10);

The main code is:

    phi=0:0.25:40*pi;
    tp=74.61;
    Ro=460;
    z=12;
    syms n m r
    b=symsum(symsum(vpaintegral(J(2*m-1,:)*380/Ro,r,[0 Ro]).*interp1(v,J(2*m-1,:),s(2*m-1,2*n-1)*380/Ro).*cos((2*n-1)*pi*phi/tp)*(exp(s(2*m-1,2*n-1)*z/Ro)+exp((2*l-z)*s(2*m-1,2*n-1)/Ro))./((2*n-1)*pi*interp1(v,J(2*m,:).^2,s(2*m-1,2*n-1))),n,1,10),m,1,10);
    B=double(b);
    plot(phi,B);

The issue is how to refer on a specific cell of the S array and the Bessel waveform array using the symbolic variables m and n. I know that with symbolic variables this cannot really happen so i am looking for a convertion command in order to keep the symbolic summation.

[Walter Roberson]

Note that https://www.mathworks.com/matlabcentral/fileexchange/48403-bessel-zero-solver is an improved version of the besselzero function you are using.

[Walter Roberson]

The integral in your equation is ambiguous. ρ is an input to the function, but the integral is . Notice the \rho in the argument to and the multiplication by ρ . It is not possible to tell at that point which of the two ρ is the variable of integration for and which is the input to the function

[Georgios Skarmoutsos]

Thanks for noting the new function for the bessel zeros. I use ρ=380 embedded on the code when I use interp1 to find specific values: ->

interp1(v,J(m,:),s(m,n)*380/Ro) % Ro=460

For the integral I use this command:

vpaintegral(interp1(v,J(m,:),s(m,n)/Ro)*r,r,[0 Ro])

Where I haven't put a value on the variable and it varies as a symbolic variable. For every different Bessel root we get another definite integral approximation. The only problem that is left on the code is how to index zeros on the formula using s(m,n) since matlab compiler doesn't allow the symbolic variable indexing. I tried using two for loops but I get funny results. I found this answered question https://uk.mathworks.com/matlabcentral/answers/300634-array-index-in-symsum?s_tid=mwa_osa_a where you also solve the same problem using sum but I don't understand how to apply the sum function to get the final f(φ) waveform. You mean something like this?

syms r
phi=0:0.2:200;
b=zeros(size(phi));
n=zeros(20,1);
m=zeros(20,1);
for i=1:2:40
    n(i)=i;
    m(i)=i;
end

b=sum(sum((8*1.37752/(n*pi*Ro^2)*(vpaintegral(interp1(v,J(m,:),380*s(m,n)/Ro)*r,r,[0 Ro])*sin(0.5*n*pi*ap)*cos(pi*n*phi/tp)*interp1(v,J(m,:),s(m,n)*380/Ro)*(exp(s(m,n)*z/Ro)+exp((2*l-z)*s(m,n)/Ro))))./((interp1(v,J(m+1,:).^2,s(m,n)*380/Ro)*(exp(s(m,n)*hm/Ro)*(1-coth(s(m,n)*hm/Ro)*1.04)+exp((2*28.5-hm)*s(m,n)/Ro)*(1+coth(s(m,n)*hm/Ro)*1.04))))));

[Walter Roberson]

What is the value of l as in exp(2*l-z) ?

Where did the coth() term come from?

[Georgios Skarmoutsos]

Sorry the full equation is this one, I neglected some constant terms. The full equation is shown below. The constants ap, hm, l, Ro, mr, are given on the code below. coth is the hyperbolic cot which is equal with coth(x)=cosh(x)/sinh(x)=(exp(x)+exp(-x))/(exp(x)-exp(-x)). The only error that is left I think is the invalid indexing for arrays J and s.

s=zeros(100,100);
for u=1:100                 % s is the array with the Bessel function zeros
    s(u,:)=bz(u,100,1);
end
v=0:0.1:9.9;
J=zeros(100,100);  %
for i = 1:100
    J(i,:) = besselj(i,Z);
end
phi=0:0.25:40*pi;
tp=74.61;
hm=10;
ap=0.75556;
l=28.5;
m0=1.25663706e-06;
mr=1.04;
p=16;
Ro=460;
R=380;
Br=1.28;
z=10+0.5*3.5;
%*interp1(v,besselj(2*m,v,s(2*m,2*n-1)*380/Ro))./((interp1(v,(besselj(2*m+1,v,s(2*m+1,2*n-1)))).^2)*(exp(s(2*m,2*n-1)*hm/Ro)+exp((2*l-hm)*s(2*m,2*n-1)/Ro))*(pi*(hm+mr*l))*)
%((exp(s(2*n-1,2*m)*z/Ro)+exp((2*l-z)*s(2*n-1,2*m)/Ro))


syms n m r
b=symsum(symsum((8*Br*vpaintegral(interp1(v,J(2*m-1,:),s(m,n)*380/Ro)*r,r,[0 Ro]).*sin((2*n-1)*pi*ap/2).*interp1(v,J(2*m-1,:),s(2*m-1,2*n-1)*380/Ro).*cos((2*n-1)*pi*phi/tp)).*(exp(s(2*m-1,2*n-1)*z/Ro)+exp((2*l-z)*s(2*m-1,2*n-1)/Ro))./((2*n-1)*interp1(v,J(2*m,:).^2,s(2*m-1,2*n-1))*pi*(mr-1)*(exp(s(2*m-1,2*n-1)))),n,1,10),m,1,10);
B=double(b);
plot(phi,B);

Answer 1

[Walter Roberson]

Symbolic variables can never be used to index an array. In particular, you are attempting J(2*m-1,:) . The ':' tells us that J is almost certainly an array, so you are trying to index using a symbolic variable.

When you have a finite set of things to index in an expression, you should not use symsum() or symprod(): instead you should create the full list of values as arrays, and sum() or prod() as appropriate.

You might need to vectorize along a 3rd or even 4th dimension to get your calculations right to end up with a 2D array as the overall result.

[Georgios Skarmoutsos]

You mean that I need to create two n and m vectors with the values of the indexes [1 3 5 7 ..] and use a double sum?

Is there any example code in mathworks which applies the sum function to plot a square wave pulse?

[Walter Roberson]

l = 4; warning('You need to set l to something appropriate!');

trunc_order = 100;
m = 1 : 2 : trunc_order;
[M,N] = ndgrid(m);
SMn = besselzero(m.', trunc_order);
SMN = SMn(:,1:2:end);

rho = 380;
z = 12;
R0 = 480;
tp=74.61;

phi = 0:0.25:40*pi;


SMN_R = SMN./R0;
part1_fun = @(r) besselj(M, SMN_R).*r;
part1 = integral(part1_fun, 0, R0, 'arrayvalued', true);
part2 = N .* pi .* besselj(M+1, SMN).^2;
part3 = besselj(M, SMN_R .* rho);
part4 = cos(N.*pi./tp.*reshape(phi,1,1,[]));
part5 = exp(SMN_R .* z) + exp((2*l - z).*SMN_R);

all_together = part1 ./ part2 .* part3 .* part4 .* part5;
f_rho_phi_z = squeeze( sum( all_together, [1 2]) );

plot(phi, f_rho_phi_z);

[Georgios Skarmoutsos]

It works perfectly!

To plot a single summation with m=n how could I modify the code?

[Walter Roberson]

If you mean one particular numeric m = n, then

squeeze(all_together(that_value, that_value, :))

If you want all of the m=n, but added together, sum of the diagonal for each theta, then perhaps the easiest approach would be

mask = repmat(diag(size(all_together,1), size(all_together,2)), 1, 1, size(all_together,3));
f_rho_phi_z = squeeze( sum( all_together .* mask, [1 2]) );

[Georgios Skarmoutsos]

I will right a summarization to make sure that I understood everything and If I am making a mistake let me know.

The sequence that the initial code computes the sum is that for m=1 adds all the terms n=1,3,5,.. then m=2 and readds the terms with index n and eventually all the m=1,3,5.. are added right?

The modification

squeeze(all_together(that_value, that_value, :))

gives a specific term of this series with m=n

The modification

mask = repmat(diag(size(all_together,1), size(all_together,2)), 1, 1, size(all_together,3));
f_rho_phi_z = squeeze( sum( all_together .* mask, [1 2]) );

computes all the terms with index m=n=1,3,5 so if I got this right it computes the following which m=n=1,3,5..:

It gives the error:

Array dimensions must match for binary array op.

Error in (line 28)
f_rho_phi_z = squeeze( sum( all_together .* mask, [1 2]) );