Implementing double series in MATLAB

[Atr cheema]

Mod note: Moved from a technical forum section, so missing the homework template
I want to write code for this double sum in MATLAB and I have written following code:

x = 100; % to calculate omega and u
l = 300; % to calculate omega
p = 10;
omegaa= x/l; deltaH = 200;
deltat = 5; % to calculate u
v = 140; % to calculate u
u = x/sqrt(v*deltat); suma = 0;
for m = 1:p;
for n = 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;
summ = suma + f;
end
end​

My problem is no matter what value of deltaH I use e1 and e2 values remain zero, may be because of p-m. How can I resolve this issue?

 

[jasonRF]

Several things:

1. e1 and e2 have no dependence on deltaH, so you shouldn't expect e1 and e2 to change with deltaH
2. You are not implementing the sum in the link. take another look and check your code. There are several errors on your part.
3. erfc(x) falls off very fast for positive x. for example, erfc(4) is about 1.5e-8, erfc(6) is about 2e-17. So if you give it large positive arguments you will get zero.

jason

 

[Atr cheema]

Yes the actual problem is the values of e1 and e2 always remaining constant.
'Summ' should be 'suma' but it won't solve problem of e1 and e2 being always zero.
Value of x is measured from field. It can not be changed to other value.

 

[Dr Transport]

First of all, you have some code errors in there, you are not calculating what you posted. secondly, the erfc function is ~10^{-100} for x > 5 so it doesn't surprise me that it come out to be 0. I inserted some parentheses and fixed it according to your image of the equation...

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;
for m = 1:p
  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)/(1*omega*sqrt(p-m));
  e2 = erfc(e2Arg);
  totalerror = e1-e2;
  f = (-1)^(n-1) * deltaH * totalerror;
  suma = suma + f;
  end
end

 

[Atr cheema]

@Dr Transport Thank you very much. One more thing. I want MATLAB to run summation for every single value of 'p', every time with a different value of 'deltaH'. It means if 'p=100' then my goal will be to calculate sum first for p=1, than for p=2 and so on till p=100, thus 100 values in total. I have added an outer for loop but final array results in an ever increasing values which should not be the case. The elements of outer array should behave somehow similar to that of deltaH. i.e bell shaped

x = 100; % to calculate omega and u
l = 300; % to calculate omega
p = 100;
omega= x/l; deltaH =sin(linspace(0.3,2.6,100));
deltat = 5; % to calculate u
v = 140; % to calculate u
u = x/sqrt(v*deltat); suma = 0;
for i=1:p
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)/(1*omega*sqrt(p-m));
  e2 = erfc(e2Arg);
  totalerror = e1-e2;
  f = (-1)^(n-1) * deltaH(i) * totalerror;
  suma(i) = suma(i) + f;
  end
end

 

[Dr Transport]

e2Arg: the denominator should have a factor of 2 not 1 and this will run faster if you use vectors not loops...

 

[Atr cheema]

@Dr Transport Having changed the denominator the result still looks unrealistic as here http://[PLAIN]https://www.dropbox.com/s/glaq5av7fzzlhaj/11.png?dl=0 . It still seems that I am adding more in every element of the output array `suma`.

 

[Dr Transport]

Look at my original reply, you are taking the difference between two extremely small numbers...

 

[Atr cheema]

@Dr Transport It seems I am not able to clear my point. You said in your initial reply that erfc is very small for x>5 but here problem is that output array "suma" is getting larger and larger which means I am adding more in every successive element of output array "suma".

 

[Dr Transport]

I'll have to interrogate it further, kinda tied up right now...

 

[Dr Transport]

I had to rework the inner loop, that is the reason your suma is increasing. I remember that your inner loop on m was not set up correctly, I fixed that in my version. when I run it on my machine, suma hits 74.59 in about 2 iterations and doesn't change.