1

We are coding with python in my numerical analysis class and I am quite new to it, so I am not very sure I did it right and would appreciate your help. With python I need to compute:

i) 0.01^2 and enter the answer rounded to 5 significant figures

ii) f(0.000122) where f(x) = cosh(x) + cos(x) - 2, and round your answer to 18 decimal places.

iii)Calculate the absolute error | f(x) - p_4(x) | at the point 0.981 if p_4(x) is the 4th-degree Taylor polynomial approximating f(x) = exp( pi*x ) with x_0 = 0

What I did:

i) My code would simply be

x=0.001**2:
    print(x)

Here I get 0.0001 but how should I round this now to 5 figures?

ii) This time I had:

import numpy as np
x = 0.000122
y = np.cosh(x) + np.cos(x) - 2
print(round(y, 18))

But the solution gives me 0.0 and this seems wrong because... again it says round it to 18 decimal places.

iii) Here I have just started with the taylor series:

import math
x = math.pi*0.981
p4 = x**0/math.factorial(0) + x**1/math.factorial(1) + x**2/math.factorial(2) + x**3/math.factorial(3) + x**4/math.factorial(4)
print(p4)

Is this ok? Where do I use that $x_0=0$? I am always unsure, when I don't use every hint. My next step then would be:

error=abs(math.e**(math.pi*0.981)-p4)

Is this ok?

Improve this question
5
  • When they say "round", they might be asking you to display that many decimal places (which is not the same thing as rounding with round()). You can round only for display purposes by using string formatting of floating point numbers. Try f"{1/3:.18f}" (with the leading f before the string - it's an f-string) in a Python version 3.6 or newer to see one third rounded to 18 decimal places. Commented Oct 28, 2019 at 13:31
  • You do not need to import numpy, the cosine in both variants is also present in the math module. // @blubberdiblub : To get the correct number of non-zero decimals, would it not be better to use the e format specifier? Commented Oct 28, 2019 at 16:26
  • @Dr.LutzLehmann that depends on what you want. As the e specifier forces e-notation (i.e. x * 10^m format), if you specify f"{0.001:.18e}", you'd get the 18 decimals after the most significant non-zero digit, whereas for f"{0.001:.18f}" you get the 18 decimals after the unit position. I assumed the latter to be the most likely interpretation, but you're of course correct that the other one could be a possible interpretation, too. Commented Nov 1, 2019 at 14:46
  • @Dr.LutzLehmann and the reason I dismissed the other interpretation is that it's not very useful to print 18 decimals (plus the decimal before the point) of a double-precision float when the decimal point is relative to the most significant digit, as IEEE double-precision is not precise enough for decimals beyond the 16th. For the f specifier it can be useful, depending on the magnitude of the value. Granted, whoever posed the question might not have known that fact. Commented Nov 1, 2019 at 15:03
  • @blubberdiblub : Yes, 18 digits is somewhat excessive, it should be enough to fix the binary pattern of the number. Without using Taylor expansions, only 8 digits are correct, so the "greater half" is garbage. Commented Nov 1, 2019 at 15:13

1 Answer 1

Reset to default
1

f(x) = cosh(x) + cos(x) - 2 has, by its Taylor expansion, a value of

[1+x^2/2+x^4/24+...] + [1-x^2/2+x^4/24+...] - 2 = x^4/12 + O(x^8)

However, the error of the evaluation of this formula at the given value x=1.22e-4 is bounded by about 4*1.1e-16. Thus the exact value is smaller than the floating point noise. Or in other words, the formula contains 2 or 3 instances of catastrophic cancellation. It is not possible to get the desired result with the required number of correct digits.

You can use the double-argument identities to alleviate this a little, 

cosh(x)-1 = 2*sinh²(x/2)
cos(x) -1 = -2*sin²(x/2)
==> f(x) = 2*sinh²(x/2) - 2*sin²(x/2)

so that

from math import sin, cos, sinh, cosh

x=1.22e-4
y = 2*(sinh(x/2)**2-sin(x/2)**2)
print("%.18e\n%.18e"%(y, x**4/12)

gives the results

1.846112113571565333e-17
1.846112133333333120e-17

With a multi-precision evaluation, one gets the more exact values for the original formula

 1.84611 21333 33333 33357 67711 91051 10045 27621 30620 39852 27994 59662 46660 18747 71781 1746e-17

and the modified formula

  1.84611 21333 33333 33357 67711 91051 10045 27621 30620 39852 27994 59662 46660 18747 71781 17468 16309 43036e-17

Even here one gets a cancellation effect in the first formula in that the second formula provides more digits.

Comparing these to the original double precision result, we see that there only 8 digits are correct in the exact formula, corresponding to the remaining catastrophic cancellation of the identical quadratic Taylor terms. The first term of the Taylor expansion has 16 correct digits, corresponding to the relative error O(x^4)~1e-16 and the floating point precision of the same size.

Improve this answer
Sign up to request clarification or add additional context in comments.

4 Comments

Thank you very much. Now your answer focuses on ii) right? Are my answers for i) and iii) correct then? And how do you now get the answer for ii) rounded by 18 decimals?
i) There is a different number of zeros in task and solution. Also, significant digits start to count only after the leading zeros. iii) That looks correct. The expansion point x_0 is the point where you evaluate the derivatives to get the Taylor coefficients and where the expansion is centered.
I really feel bad for asking but... what do you mean with "There is a different number of zeros in task and solution"? So my answer is correct right?
In the task description, you used 0.01=1e-2. In the python code you wrote 0.001=1e-3. In the result you wrote 0.0001=1e-4, which is the square of the first one, so probably you did it right in the original code.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.