MATLAB: How to sum up large values with high precision

MATLABprecisonsymbolic

Hello all together,
I would like to calculate a formula which consists of the alternating (positive/negative) summation of very large values. For small values of N (for example N=3) ist works well. Also for larger values of N, if the gaps between the "Mueh" values are large enough. For the planned application of the formula, however, the values must be quite close together and N should nevertheless be able to reach at least values of 20 or 30. Here´s my current version using symbolic functions and vpa:
% Definition of values
N=11;
Mueh=[0.83 0.84 0.85 0.86 0.87 0.89 0.90 0.91 0.92 0.93 0.94];
Muehc=num2cell(Mueh);
T=25;
% Creating symbolic variables
syms C zaehler nenner CfF MfF;
mueh=sym('m',[1,N]);
syms t;
% Symbolic computing of the C-Values
for k=1:1:N
    C(1,k)=1;
end
for i=2:1:N
    for k=1:1:N
        zaehler=1;
        nenner=1;
        for j=0:1:(i-2)
            zaehler=vpa((zaehler*mueh(N-j)),50);
        end
        for l=0:1:(i-1)
            if (l~=(k-1))
                nenner=vpa((nenner*(mueh(N-l)-mueh(N-k+1))),50);
            end
        end
        C(i,k)=vpa((zaehler/nenner),50);
    end
end
% Symbolic computing of the function (here only for n=1)
for n=1:1:1
      for j=1:1:N-n+1
          CfF(n,j)=vpa(C((N+1-n),j),50);
      end
      for l=1:1:N-n+1
          MfF(n,l)= vpa(mueh(N-l+1),50);
      end
    f(t,mueh)=sum(vpa((vpa(CfF(n,:),50)) .* vpa(exp(-vpa(MfF(n,:),50)*t),50),50))
    h(t)=f(t,Muehc{:}) 
  end
% Plot h
fplot(h, [0 5]);
The current result – although using sym and vpa looks like this:
I also tryed XSum without success.
Do you have an idea how to fix it?
Many thanks and best regards, Johannes

Best Answer

There are a LOT of problems in this code, if you are trying to obtain high precision in the result.
To start:
Mueh=[0.83 0.84 0.85 0.86 0.87 0.89 0.90 0.91 0.92 0.93 0.94];
You may THINK you have created those values exactly. But you did not. In fact, they are approximate, correct to roughly 16 decimal digits.
sprintf('%0.55f',Mueh(1))
ans =
    '0.8299999999999999600319711134943645447492599487304687500'
As you can see, what MATLAB stored was only an approximation. The same is true for the rest of that vector. You can do further computations with those numbers, but if you start with a double precision approximation, you will never have more than the 16 significant digits you stated with. A 50 digit computation is just then a waste of time and CPU cycles.
Next, I see you doing things like this:
% Creating symbolic variables
syms C zaehler nenner CfF MfF;
But then you assign them with values as
zaehler=1;
WRONG. You are thinking as if that variable is forever after defined as symbolic! Not true. Using syms is NOT a permanent declaration of class, something that sticks forever for that variable.
clear
syms X
whos X
Name      Size            Bytes  Class    Attributes
X         1x1                 8  sym
So MATLAB knows that X is symbolic. However, if I now assign a double into X,
X = 3;
whos X
Name      Size            Bytes  Class     Attributes
X         1x1                 8  double
X is now a double precision number.
When you make an assignment like that, X gets overwritten. Everything about X gets overwritten.
I also see you trying to index into an array using a symbolic index, t. You can't do that.
Let me next see if I can re-write your code so it will be more accurate...
I'll use my own hpf tools.
% Definition of values
DefaultNumberOfDigits 100 5
N=11;
mueh=hpf(83:94,100)/100; % These are now exact.
% Compute C in HPF form
C = zeros(N,N,'hpf');
C(1,:) = 1; % small integers are exact, as flints
for i=2:1:N
  for k=1:1:N
      zaehler=hpf(1);
      nenner=hpf(1);
      for j=0:1:(i-2)
          zaehler=zaehler*mueh(N-j);
      end
      for l=0:1:(i-1)
          if (l~=(k-1))
              nenner=nenner*(mueh(N-l)-mueh(N-k+1));
          end
      end
      C(i,k)=zaehler/nenner;
  end
end
So the above computation of C is done correctly, carrying 100 digits, with 5 guard digits as insurance. That is what the statement
DefaultNumberOfDigits 100 5
does for you. The 5 is the number of extra guard digits, so you could specify more if that would be useful.
You should be able to do the remainder.
Related Question