Finding erf using Maclaurin series breaks down around x = 5

After I got around 2.9850991e+40 instead of 1 at x=10, I thought it may be related to precision, so I tried this code:

X = [0.5, 1, 5, 10]
for i = 1:length(X)
    x = X(i)
    fprintf('Actual: %#.8g', erf(x))
    ErrorFunctionNew(x)
end

function ErrorFunctionNew(x)
SUM = vpa(0,1000);
for n = 0:100
    N(1,n+1) = vpa(sym((((-1)^n)*(x^(2*n+1)))/(factorial(n)*(2*n+1))),1000);
end
%disp(N)
SUM = vpa(sum(N,'all'),1000);
base = vpa(sym(2/sqrt(pi)),1000);
Answer = base * SUM;
fprintf('Erf: %#.8g', Answer)
end

But it did not really help, what am I doing wrong, or is it just not possible with this formula?

Here is the source:

Update: The professor sent this code back after I messaged him for help with my revised code. (Question was homework, and deadline is still weeks away, but I guess he was leniant on this one)

x = 0.5;
Erf = 0;
sign = 1;
fac = 1;
for n = 0:500
    Erf = Erf + sign*(x^(2*n+1))/(fac*(2*n+1));
    sign = -sign;
    fac = fac*(n+1)
end
Erf = Erf*2/sqrt(pi)

[Jan]

The faculty, power and exp are very expensive operations. Avoiding (-1)^n by a repeated multiplication by -1 is cheap. This works for x^(2*n+1) also:

x   = 0.5;
Erf = 0;
sgn = 1;
fac = 1;
xp  = x;
xpm = x * x;
for n = 0:500
    % Sum of: (-1)^n * x^(2*n+1) / (factorial(n) * (2*n+1))
    Erf = Erf + sgn * xp / (fac * (2 * n + 1));
    sgn = -sgn;             % (-1)^n
    xp  = xp * xpm;         % x^(2n+1)
    fac = fac * (n + 1);    % faculty(n)
end
Erf = Erf * 2 / sqrt(pi);
Erf - erf(x)

ans = 1.1102e-16

The sum until 400 is a waste of time: fac is Inf for n=170 already and the following terms of the sum are 0, so insert:

if isinf(fac), break; end

For x=7 the method fails due to overflow.

[Phuwanut Pataratawinun]

That was what I thought too, and I messaged them, but they kind of avoided answering the question at all costs so I guess I will just use the method with vpa.

Thank you for your help though.

Answer 1

[Jan]

Your code is almost correct and working, except for: factorial(n). When n is a double, this tends to overflow soon. So calculate the factorial symbolically also:

x = sym(7);
digits(200) 
for n = 0:200
    sn       = sym(n);
    N(1,n+1) = vpa(sym((((-1)^n) * (x^(2*n+1))) / (factorial(sn) * (2*n+1))));
    %                                                        ^^
end
SUM    = vpa(sum(N));
base   = vpa(sym(2/sqrt(pi)));
Answer = vpa(base * SUM)

Answer = 

0.99999999999999998640926093836298196491479643119628317876340337678863707493643361627341573937940979508185850406740520056891281700621873271964658015774990123398226165203422093626921128793078138097349342

double(Answer - erf(x))

ans = -1.3591e-17

Now increasing the number of digits or the length of the sum does not increase the accuracy of the output. But why? Again: sym(2 / sqrt(pi)) is calculated with doubles at first.

base   = vpa(sym(2) / sqrt(sym(pi)));
Answer = vpa(base * SUM)

Answer = 

0.99999999999999999999995816174392220585654792132738314163365262144769685562329269344004967730248854760989606890533717031954476805761887825825106673628695780003101195720085909248130607112982333829892876

double(Answer - erf(x))

ans = 5.3406e-40

Now digits(1000) and for n=0:400 increase the accuracy such that the deviation is 2e-196.

[Torsten]

What about the recursion

N(1,1) = x;
N(1,n+1) = N(1,n) * (-1) * x^2 * (2*n-1)/(2*n+1) * 1/n ;

?

No factorial needed.

Answer 2

[David Hill]

ERF=@(x)round(2/sqrt(pi)*integral(@(t)exp(-t.^2),0,x),8);%much easier using the integral definition

[David Hill]

Alternatively (although it runs into floating point problems),

function e = ERF(x)
e=0;
for n=0:1000
    e=e+x/(2*n+1)*prod(-x^2./(1:n));%after modifying the equation slightly (see https://en.wikipedia.org/wiki/Error_function) 
end
e=2*e/sqrt(pi);
end

[Torsten]

For x=10, it gives the same problem as in the OP's code.

[David Hill]

Yes, integral solution works.

Answer 3

[John D'Errico]

Do you expect that all infinite series are convergent? And sometimes, even a series that is theoretically convergent for all x, do not work, because of numerical problems, which would force you in theory to work in extremely high precision.

https://en.wikipedia.org/wiki/Error_function

Note the comment in there about how the erf series is famously known for poor convergence even for x greater than only 1. A problem is the denominator coefficients do not go to zero quickly enough to kill off the numerator, until the number of terms grows very large. And since the terms have alternating signs, that kills you with massive subtractive cancellation. (If you don't know what that expression means, do some reading online.)

The problem is not in your code. It is the series itself. Having said that, are there tricks one can use? Well, yes. There are other series that are not technically Taylor series, but you are asking how to deal with erf in context of the Taylor series. As I recall, you can also work with the complimentary error function for large values of x. I'm sure there are other tricks one can do, but your question is purely in terms of the Taylor series.

Sorry, but you were effectively given this assignment to learn that things can get nasty, not to actually get a good answer.

Look at the terms, even for small x.

x = sym(5);
n = (0:10)';
vpa((-1).^n.*x.^(2*n+1)./(factorial(n).*(2*n+1)))

ans = 

Do you see that each term is growing in size quickly? Eventually, the x^(2*n+1) will dominate. For example....

N = [20 40 60 80 100]';
vpa(x.^(2*N+1)./(factorial(N).*(2*N+1)))

ans = 

So finally, after 100 terms, they are starting to get small. But if you look at the intermediate terms, you will see that massive subtractive cancellation kills your chances. If x is larger yet, like 10, things go to hell way worse. Again, this is a normal problem, given to students purely to see them squirm when they think they got things wrong. I can think of a few others that are classically nasty. Thank your instructor.

Seriously, if you really wanted to evaluate the error function, there are better approximation forms than the series you were told to use here. I would start by loooking in say Abramowitz & Stegun, where you might find perhaps a Pade approximant that will do better. And then I would look at the classic by JF Hart, "Computer Approximations", which surely has many such approximations.

[Phuwanut Pataratawinun]

First of all, thank you for the quick and detailed answer.

I know there are many major problems, but if I just needed the final answer I could just get it from erf(x). Though that is beside the point. Is there a conceivable way to make code that computes erf based on the given formula from the problem(see attached picture in the original question) in MATLAB?

Here is my revised testing code:

clearvars
x = sym(7)
digits(1000) 
fprintf('Actual:    %#.8g', erf(x))
for n = 0:200 
    N(1,n+1) = vpa(sym((((-1)^n)*(x^(2*n+1)))/(factorial(n)*(2*n+1))));
end
M = zeros(length(N),1);
for i = 1:length(N) % Convert to double for easier viewing 
    M(i) = N(i);    % with variables viewer to see if n is big enough
end                 % or if digits() needs an increase
SUM = vpa(sum(N,'all'));
base = vpa(sym(2/sqrt(pi)));
Answer = vpa(base * SUM)
fprintf('Erf Func: %#.8g', Answer)

Looking through M with the variable viewer I see the first flat 0 with no decimals is at position 171, and viewing N(1,171) with the variable viewer confirms this, so I think n of 200 should work for x = 7.

Also, I see the biggest number is at position 48, and viewing N(1,48) with the variable viewer gives -7832901140721959410.142200 ... so I think digits(1000) should work also.

Problem is erf(7) gives 1.0000000 which I should get close to, but instead I get 2972.86239339 ... .

Where is the inaccuracy(ies) coming from or where did I go wrong?

[Torsten]

I think John D'Errico is right:

This is a problem you should learn from, not solve it.

[Jan]

@John D'Errico: As soon as the factorial is evaluated with a high precision also, the strange result vanishes and the sum converges:

factorial(n)  -->  factorial(sym(n))