Normalized to unity by summation

I have the following code where I am trying to get unity after performing a summation. But it is not giving me unity using '' symsum''. Is there any alteranative way to execute the summation with a complicated expression. Pl somebody help me to solve that. Here my code is given below. Although theoretically it is giving ''one'' for n ranging from 1 to inf.

clc;
syms  n 
x=1.0;
y=0.5;
g=0.2;
l=0.1;


om=sqrt(((x).^2)-(4.*(g.^2)));

mu=sqrt((x+om)./(2.*om));

nu=((x-om)./(2.*g)).*mu;

eta=(((l)./((2.*g)+x)).*(1+((x-om)./(2.*g)))).*mu;


%% n dependency start

En=((n+(1./2)).*om)-(x./2)-(((l).^2)./((2.*g)+x));

Enn=((n-(1./2)).*om)-(x./2)-(((l).^2)./((2.*g)+x));

Dn=(y./2).*(exp(-2.*((eta).^(2)))).*(laguerreL(n,(4.*((eta).^2))));     

Dnn=(y./2).*(exp(-2.*((eta).^(2)))).*(laguerreL((n-1),(4.*((eta).^2)))); 

Em=En - Dn;

Ep=Enn + Dnn;

eps=(Ep-Em)./2; 

Deln=(eta.*y./sqrt(n)).*exp(-2.*(eta.^2)).*laguerreL((n-1),1,(4*(eta^2)));

xn=sqrt(((eps).^2)+((Deln).^(2))); 

zetap=sqrt(((xn.^2)+(eps.^2))./(2.*xn));

zetam=sqrt(((xn.^2)-(eps.^2))./(2.*xn));

z= 1i.*(mu-nu).*eta./sqrt(2.*mu.*nu);

a1=(zetap./sqrt(factorial(n-1))).*((- nu./(2.*mu)).^(-1./2)).* hermiteH(n-1, z);

b1=(Deln./abs(Deln)).*(zetam./sqrt(factorial(n))).* hermiteH(n, z);

a2=(zetam./sqrt(factorial(n-1))).*((-nu./(2.*mu)).^(-1./2)).* hermiteH(n-1, z);

b2= (Deln./abs(Deln)).*(zetap./sqrt(factorial(n))).*hermiteH(n, z);


c0= -(1./sqrt(2.*mu)).*exp(-((eta.^2)./2)+ ((nu.*(eta).^2)./(2.*mu)));

cp= -c0.*((-nu./(2.*mu)).^(n./2)).*(a1 - b1);

cm= -c0.*((-nu./(2.*mu)).^(n./2)).*(a2 + b2);

c=((abs(cp)).^2) + ((abs(cm)).^2);

c1=(abs(c0)).^2;

out= ((abs(c0)).^2) + symsum(c,n,1,100);

vpa(out,5)


%    Actually I need to get    c1 + sum_{n=1}^{inf} c = 1. Is it possible to get by using symsum??

Answer 1

[Walter Roberson]

You can use

out= ((abs(c0)).^2) + symsum(c,n,1,inf);

and that will be faster than for upper bound 100. The symsum() will stay unevaluated, and will be converted into numeric form when you vpa().

Alternately you could use

out= ((abs(c0)).^2) + sum(subs(c,n,1:100));

This will be a bit slow. I would suggest instead

guard_digits = 10;
out= ((abs(c0)).^2) + sum(vpa(subs(c,n,1:100), guard_digits);

Any of these methods will give the same result to within round-off error -- and none of them will give 1.0 for the final result, not even the one with symsum() to infinity. I recommend that you re-check the formulaes.

[AVM]

@walter: Thanks for your help. Yes, I did a mistake in one of the expression. Now it exactly 'one'