implementing double $\sum_{m=1}^{p}\sum_{n=1}^{\infty }$ sum in MATLAB

Question

I want to calculate following double sum: $$\sum_{m=1}^{p}\sum_{n=1}^{\infty }(-1)^{n-1}\Delta H.{erfc\frac{2(n-1)u+u\gamma}{2\gamma \sqrt{p-m}}+erfc\frac{2nu-u\gamma}{2\gamma \sqrt{p-m}} }$$ I have written following code in MATLAB

x = 45;         % to calculate omega and u  
l = 300;          % to calculate omega
p = 100;  
omega= x/l;   deltaHvec = sin(linspace(0.3,2.6,100));
deltat = 5;      % to calculate u
v = 40;          % to calculate u
u = x/sqrt(v*deltat); suma = 0; 
for i = 1:p
    deltaH=deltaHvec(i);
    for m = 1:i
        for n = 1:1000
            e1Arg=(2*(n-1)*u + u*omega) / (2*omega*sqrt(p-m));
            e1=erfc(e1Arg);
            e2Arg = (2*n*u - u*omega)/(2*omega*sqrt(p-m));
            e2 = erfc(e2Arg);
            totalerror = e1+e2;
            f= (-1)^(n-1) * deltaH * totalerror;

            summ(i) = summ(i) + f;
        end
    end

end

summ;
plot(1:100,summ); hold on; plot(1:100,deltaHvec,'k');

I want to have 100 values of summ i.e One for every different value of 'p' every time using one specific value of 'deltaHvec'. I mean for first first value of 'summ', p=1 and deltaH=deltaHvec(1). For second value of 'summ' 'p=2' and deltaH=deltaHvec(2). The problem is final graph of 'summ' should be equal or less than amplitude of 'deltaHvec'. But I am getting quite higher values of 'summ' which represents that I am adding extra values in 'summ' with every successive step.

Answer 1

If I understood your problem correctly, you just want to avoid the cases where p=m in your double loops. You can test that with a "if" function but using parallelization is better.

x = 100;         % to calculate omega and u  
l = 300;          % to calculate omega
p = 10;  
omega = x/l;
deltaH = 200;
deltat = 5;      % to calculate u
v = 140;          % to calculate u
u = x/sqrt(v*deltat); suma = 0; 

[m,n] = ndgrid(1:p,1:1000);

e1 = erfc(2*(n-1)*u + u*omega ./ (2*omega*sqrt(p-m)));
e2 = erfc(2*n*u - u*omega./(1*omega*sqrt(p-m)));
totalerror = e1.*e2;
f = (-1).^n-1 * deltaH .* totalerror;

f(eye(size(m))==1) = 0;
f = sum(sum(f,1),2);

I've not tried to run it since the definition of some variables are missing.