problem with vpa and symsum

Hi i make the fumtiom

function z= Escattheory(r,phi)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
syms  kk
[f,N,Nc,a,ra,k0,Z0] = parameter();
E0=1;
Z1=k0.*r;
Z2=k0.*ra;
factor1=besselj(kk,Z2).*besselh(kk,2,Z1)./besselh(kk,2,Z2)
factor2=(-1j).^kk.*cos(kk.*phi)
factor=-E0.*factor1.*factor2.*en(kk);
%z=vpa(symsum(factor,k,0,inf),3)
f_all=symsum(factor,kk,0,300);
y=vpa(f_all,3)
end

when i call the function Escattheory(r,phi) with argument r=ra=1 and phi =0 the resultas is - 1.22 + 5.05e-5i

if i use the sum like this

function z= Escattheory_test(r,phi)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
Escattheory_test(r,phi)
    function   y=Escattheory_test(r,phi)
        %UNTITLED Summary of this function goes here
        %   Detailed explanation goes here
        [f,N,Nc,a,ra,k0,Z0] = parameter();
        E0=1;
        Z1=k0.*r;
        Z2=k0.*ra;
        kk=0:200
        factor1=besselj(kk,Z2).*besselh(kk,2,Z1)./besselh(kk,2,Z2)
        factor2=(-1j).^kk.*cos(kk.*phi)
        factor=-E0.*factor1.*factor2.*e_n(kk);
        %z=vpa(symsum(factor,k,0,inf),3)
        y=sum(factor)
    end

the reuslts is -1-j0 what happen ??

thank you

[Torsten]

You forgot to include the parameter function and the e_n (en) function.

[george veropoulos]

thank the parameter function is

function [f,N,Nc,a,ra,k0,Z0] = parameter()
%UNTITLED Summary of this function goes here
c0=3e8;
Z0=120.*pi;
ra=1;
a=0.92;
N=30;
Nc=30;
f=300e6;
lambda=c0./f;
k0=2*pi./lambda;
end

the function e_n is

 
function y=e_n(k)
if  k==0
    y=1
elseif k~=0
y=2
end
end

give only numerical value

Answer 1

[Torsten]

The numerical bessel functions are not able to compute your sum correctly (for higher values of k, they return NaN) (see below). That's the reason why you get different results.

r = 1;
phi = 0;
Escattheory(r,phi)
Escattheory_test(r,phi)

factor = 
  Columns 1 through 9

  -0.4406 - 0.0000i   0.0000 - 0.4248i  -0.5758 + 0.0000i   0.0000 - 0.0582i  -0.6314 - 0.0000i   0.0000 + 0.7456i   0.5554 + 0.0000i   0.0000 - 0.3150i  -0.1466 + 0.0000i

  Columns 10 through 18

   0.0000 + 0.0582i   0.0203 - 0.0000i   0.0000 - 0.0063i  -0.0018 + 0.0000i  -0.0000 + 0.0005i   0.0001 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 19 through 27

   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 28 through 36

   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i

  Columns 37 through 45

  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i

  Columns 46 through 54

   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 55 through 63

   0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 64 through 72

   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i

  Columns 73 through 81

  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i

  Columns 82 through 90

  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i

  Columns 91 through 99

   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 100 through 108

   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i

  Columns 109 through 117

  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i

  Columns 118 through 126

   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 127 through 135

   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 136 through 144

   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i

  Columns 145 through 153

  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i

  Columns 154 through 162

   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 163 through 171

   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 172 through 180

   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i

  Columns 181 through 189

  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i

  Columns 190 through 198

   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 199 through 207

   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 208 through 216

   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i

  Columns 217 through 225

   0.0000 + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i

  Columns 226 through 234

   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi

  Columns 235 through 243

      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i

  Columns 244 through 252

   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi

  Columns 253 through 261

      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i

  Columns 262 through 270

   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi

  Columns 271 through 279

      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i

  Columns 280 through 288

   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi

  Columns 289 through 297

      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i

  Columns 298 through 301

   0.0000 +    NaNi      NaN + 0.0000i   0.0000 +    NaNi      NaN + 0.0000i

ans = NaN +    NaNi

function z= Escattheory(r,phi)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
syms  kk
[f,N,Nc,a,ra,k0,Z0] = parameter();
E0=1;
Z1=k0.*r;
Z2=k0.*ra;
factor1=besselj(kk,Z2).*besselh(kk,2,Z1)./besselh(kk,2,Z2);
factor2=(-1j).^kk.*cos(kk.*phi);
factor=-E0.*factor1.*factor2.*e_n(kk)
%z=vpa(symsum(factor,k,0,inf),3)
f_all=symsum(factor,kk,0,300);
z=vpa(f_all,3);
end

function   y=Escattheory_test(r,phi)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
[f,N,Nc,a,ra,k0,Z0] = parameter();
E0=1;
Z1=k0.*r;
Z2=k0.*ra;
kk=0:300;
factor1=besselj(kk,Z2).*besselh(kk,2,Z1)./besselh(kk,2,Z2);
factor2=(-1j).^kk.*cos(kk.*phi);
factor=-E0.*factor1.*factor2.*en(kk)
%z=vpa(symsum(factor,k,0,inf),3)
y=sum(factor);
end

function [f,N,Nc,a,ra,k0,Z0] = parameter()
%UNTITLED Summary of this function goes here
c0=3e8;
Z0=120.*pi;
ra=1;
a=0.92;
N=30;
Nc=30;
f=300e6;
lambda=c0./f;
k0=2*pi./lambda;
end
function y=en(k)
if  k==0
    y=1;
else
y=2;
end
end

function y=e_n(k)
y = piecewise(k==0,1,k~=0,2);
end

[george veropoulos]

thank you for your answer ! What is the best solution to keep a suficient number of terms where the bessel values are finite ? or i find other matlab function for besselj ?

When i keep 200 term in sum or for loop the reults is -1+j0 that is correct

[Torsten]

I don't know - I'd use the symbolic solution in this case:

f_all=symsum(factor,kk,0,Inf);
z=vpa(f_all);

[george veropoulos]

According. To. Physics of the problem the correct value is. From for loop or sum -1+j0

[Torsten]

According to MATLAB, the result is -1.

[george veropoulos]

The reults -1 prodeuce from the for loop or sum not from the syssum

[Torsten]

The reults -1 prodeuce from the for loop or sum not from the syssum

No.

r = 1;
phi = 0;
format long
Escattheory(r,phi)

ans = 
 -1.000000000000000 + 0.000000000000000i

function z= Escattheory(r,phi)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
syms  kk
[f,N,Nc,a,ra,k0,Z0] = parameter();
E0=1;
Z1=k0.*r;
Z2=k0.*ra;
factor1=besselj(kk,Z2).*besselh(kk,2,Z1)./besselh(kk,2,Z2);
factor2=(-1j).^kk.*cos(kk.*phi);
factor=-E0.*factor1.*factor2.*e_n(kk);
%z=vpa(symsum(factor,k,0,inf),3)
f_all=symsum(factor,kk,0,Inf);
z=double(f_all);
end

function [f,N,Nc,a,ra,k0,Z0] = parameter()
%UNTITLED Summary of this function goes here
c0=3e8;
Z0=120.*pi;
ra=1;
a=0.92;
N=30;
Nc=30;
f=300e6;
lambda=c0./f;
k0=2*pi./lambda;
end

function y=e_n(k)
y = piecewise(k==0,1,k~=0,2);
end

[george veropoulos]

i my matlab the reuslats form syssum is

- 1.2202769085399344622768816507208 - 7.6237530221056436364566359163767e-36i

>>

[george veropoulos]

with en is define from picewish the redults is -1+j0

[Torsten]

With

function y=e_n(k)
if  k==0
    y=1;
else
y=2;
end
end

MATLAB evaluates the sum

-2*sum_{k=0}^{k=Inf}(-i)^k J_k(2*pi)

with

function y=e_n(k)
y = piecewise(k==0,1,k~=0,2);
end

MATLAB evaluates the sum

-J_0(2*pi) -2*sum_{k=1}^{k=Inf}(-i)^k J_k(2*pi)

I can't exactly explain why. Maybe if-statements for the coefficients of a symbolic sum can only be implemented using a symbolic function (piecewise in this case).

[george veropoulos]

yes i think the piesewish is symbolic maybe working better in symsum

[george veropoulos]

in the for loop i use the if to define the e_n(k)

[Torsten]

As said, I don't know precisely why "piecewise" works with the symbolic solution while your if-construct gives a wrong answer.

My guess is that your function "e_n" is called with a symbolic (unspecified) variable (kk) that is not identified as "0" (which would give 1 for y), but as "else" (which always gives 2 for y).

But you might want to contact MATLAB support for a definite answer:

https://uk.mathworks.com/support/contact_us.html

[george veropoulos]

thank! i sent a report about this in matlab