MATLAB: Implementing a sum with summands of unequal spacing inside an integral

fourier seriesMATLABsum inside integralsum of unequal spacing

Dear all,
for a university study, I need to implement the following equations which gives me some headaches as I'm not a MATLAB expert.
For reasons of good scientific practice, it shall be noted that the following equations are citated from: "Analytical calculation of D- and Q-axis Inductance for interior Permanent Magnet Motors Based on Winding function Theory", Liang et al. in Energies 2016, 9, pp.580-591.
where:
N, k_wv, p, I, Rs, Rr and L, are constants. Same holds for U_d1, U_d2, alpha, beta and mu_0 (see below). Besides, B_d is a combination of multiple Fourier-series:
There are three big issues that I'm struggeling with:
1) While j=1,…,n, for v holds v=1,5,7,11. Therefore I wasn't able to use the "symsum" command as the definition of lower and upper bound for v would include that every integer within 1-11 would be considered.
2) I'm not sure if I can just say (1:j) meaning j sould go from 1 to J and define J=100 later on because…
3) The error message I'm receiving is " Index exceeds number of array elements (1)" in the line "DataStorage(1,i)=B_d…
Although I'm already searching for a while to solve 3), I'm almost sure that my approach isn't purposeful. So any help is highly appreciated!
Finally, please notice that I also would highly appreciate every piece of advice considering only a part of this mess here (e.g. concerning the "v-problem"). Thanks in advance.
This is what I tried:
%% Parameter
p=4;N=32;I=1;alpha=0.742;beta=0.886;Rs=0.091;Rr=0.090;mu_0=4*pi*10^(-7);
L=0.15;k_w=0.95;W=[1,5,7,11];U_d1=6.0368;U_d2=3.3872;
%% Calculation
syms theta j
B_d=@(v,j)(3*N*k_w*I/(v*p*pi)*cos(v*theta)-U_d1*4/pi*sum(sin((1:j)*alpha*pi/2)/(1:j)*cos((1:j)*theta))-U_d2*4/pi*sum((sin((1:j)*pi*beta/2)-sin((1:j)*alpha*pi/2))/(1:j)*cos((1:j)*theta)))*mu_0/(Rs-Rr);
summand=@(v,theta)sum(2*N*k_w/(p*pi)*cos(v.*theta)*B_d*(Rs+Rr)/2*L);
for i=1:length(W) %Solve the sum over v and j in the expression for B_d; that should bring a function of theta only.
    v=W(i);
    DataStorage(1,i)=B_d(v(i),100);
end
B_d=sum(DataStorage, 'all'); %Does this work even there's still a variable in B_d (theta)?
for i=1:length(W)
    v=W(i);
    L_md=integral(summand(v,theta),0,2*pi)./I;
end

Best Answer

Hi Richard,
I don't claim to be a professional although I did use Matlab at work, and here is how I would do this. The main thing has nothing to do with Matlab. It's just to recognize that the integral is from 0 to 2 pi and contains cos(nu*theta) in the integrand and another cos(j*theta) or cos(nu*theta) [separate sum so this is an independent nu from the first one] in Bd. Therefore orthogonality applies, so all that survives are the factors of cos(nu*theta) in the integrand times the factors of cos(nu*theta) in Bd, both for the same nu. Of the sum j=1 to 100, only j = [1 5 7 11] matter.
Integration 0 to 2pi gives an extra factor of pi since the average value of cos^2 = 1/2. The following result is done numerically and agrees with yours (caveat below).
p=4;N=32;I=1;alpha=0.742;beta=0.886;Rs=0.091;Rr=0.090;mu_0=4.*pi.*10^(-7);
L=0.15;k_w=0.95;W=[1,5,7,11];U_d1=6.0368;U_d2=3.3872;
nu = W;   % notation
% include everything except cos(nu*theta) or cos(j*theta) factors
intgr = (2*N*(k_w/(p.*pi))*((Rr+Rs)/2)*L/I)./nu;     
% toss in overall constant for Bd later on
Bd  =  3*N*(k_w/(p.*pi))*I./nu  -(4*U_d1/pi)*sin(nu*alpha*pi/2)./nu ...
       -(4*U_d2/pi)*(sin(nu*beta*pi/2)-sin(nu*alpha*pi/2))./nu;
result = pi*sum(intgr.*Bd*(mu_0/(Rs-Rr)))
The caveat is, you seem to be missing a factor of 1/nu in the integrand expression. Insert that, and the methods agree.
Related Question