Taylor's Approximation playing up

I'm trying to approximate a function (e^x) to a 10th order approximate about x = 0. I have made my code compatible with anonymous functions and it works for the most part. When I approximate e^x to a 8th order, it gives the correct answer, however when I go higher than the 8th order, the answer gets weird.

Here's the code:

    Func = @(x) exp(x);
    a = 0;
    N = 10;

    FuncApprox = 0;
    for i = 0:N
        syms x
        f_derrived = matlabFunction( diff(Func(x),i) );
        FuncApprox = FuncApprox + ( f_derrived(a)/factorial(i) )*( x-a )^i;
    end
    disp(FuncApprox)

When I run the code, this is what I get:

(1301357606610903*x^10)/4722366482869645213696 + (1626697008263629*x^9)/590295810358705651712 + x^8/40320 + x^7/5040 + x^6/720 + x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1

When I should get:

x^10/3628800 + x^9/362880 + x^8/40320 + x^7/5040 + x^6/720 + x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1

Answer 1

[Jan]

The result is mathematically correct, but the display could be simplified. 1301357606610903 / 4722366482869645213696 equals 1 / 3628800. Now it matters, what exactly your question is. Do you want to know how to simplify the output or why the values are displayed in this way?

Maybe simplifyFraction helps. ([EDITED] Nope, it does not.)

Move the factorial out of the symbolic part to the numerical part:

FuncApprox = FuncApprox + f_derrived(a) * (x - a)^i / factorial(i)

This works also, when i! exceeds flintmax.

[EDITED] Another solution:

FuncApprox = FuncApprox + (f_derrived(a) / factorial(sym(i))) * (x - a)^i;

[Ameer Hamza]

"1301357606610903 / 4722366482869645213696 equals 1 / 3628800"

there is a small difference caused by calculation in floating points

digits(20)
vpa(str2sym('4722366482869645213696/1301357606610903'))
ans = 
  3628800.000000000313

also, a numerator ending in 3 and denominator ending in 6 can never simply to a fraction like that.

[Jan]

Ameer, you are right. 3628800 cannot be 4722366482869645213696 / 1301357606610903 obviously. I was tired.

Answer 2

[Ameer Hamza]

You are getting this error because floating point precision issues. The result of (f_derrived(a)/factorial(i)) is initially a double value, which later becomes a sym after multiplying with ( x-a )^i. To prevent this calculation from happening in floating point numbers, convert the coefficient to sym explicitly like this (f_derrived(a)/sym(factorial(i))). This will prevent the floating point calculation and hence loss of precision. The following will give correct answer

Func = @(x) exp(x);
a = 0;
N = 10;
FuncApprox = 0;
for i = 0:N
  syms x
  f_derrived = matlabFunction( diff(Func(x),i) );
  FuncApprox = FuncApprox + ( f_derrived(a)/sym(factorial(i)) )*( x-a )^i;
end
disp(FuncApprox)

Result:

x^10/3628800 + x^9/362880 + x^8/40320 + x^7/5040 + x^6/720 + x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1

For more details of why you lose precision, refer to the answer on this question: https://www.mathworks.com/matlabcentral/answers/396433-i-want-to-do-this-2-38987696573583658648235686-5-6387649837698237687649879932-and-i-need-all-digits