by: Tomasz Romańczukiewicz
web: th.if.uj.edu.pl/~trom/
rm B-2-03
Most important problems that numerical methods have to deal with are:
Note that euclidian norm can be used instead, which is smaller but more expensive
c = np.float16(2.005)+np.float16(2.004) # ~ 4 digit accuracy
exact = 4.009
error = abs(c-exact)
print ("c={:.4e}, error={:.4e}, relative err={:.3f}%".format(c, error, error/c*100))
c=4.0078e+00, error=1.1875e-03, relative err=0.030%
c = np.float16(2.005)-np.float16(2.004) # ~ 4 digit accuracy
exact = 0.001
error = abs(c-exact)
print ("c={:.5e}, error={:.5e}, relative err={:.3f}%".format(c, error, error/c*100))
c=1.95312e-03, error=9.53125e-04, relative err=48.800%
Note 49% relative error. Only the first digit is certain!
Subtracting close numbers is the worst numerical operation
Note: Errors are added but values are subtracted. This can lead to a huge loss of accuracy (many orders of amplitude)!
c = np.float16(2.005)*np.float16(2.004) # ~ 4 digit accuracy
exact = np.float64(2.005)*np.float64(2.004)
error = abs(c-exact)
print ("c={:.5e}, error={:.5e}, relative err={:.3f}%".format(c, error, error/c*100))
c=4.01953e+00, error=1.51125e-03, relative err=0.038%
Recall the example of quadratic equation
a = 0.2; b=15; c=0.2
Δ = b**2-4*a*c
x1 = x1_exact = (-b+np.sqrt(Δ))/(2*a)
x2 = x2_exact = (-b-np.sqrt(Δ))/(2*a)
print ("x1 = ", x1, "x2 =", x2)
Δ = np.float16(b**2-4*a*c)
x1 = (-b+np.sqrt(Δ))/(2*a)
x2 = (-b-np.sqrt(Δ))/(2*a)
print ("x1 = ", x1, "x2 =", x2)
print ("Relative errors: Ex_1={:.3f}%, Ex_2={:.3f}%".\
format(100*(x1-x1_exact)/x1_exact, 100*(x2-x2_exact)/x2))
x1 = -0.01333570454687738 x2 = -74.98666429545312 x1 = -0.01953125 x2 = -74.98046875 Relative errors: Ex_1=46.458%, Ex_2=-0.008%
$x_1$ has a 56% error when accuracy of $10^{-4}$ is used wheras $x_2$ only 0.08%. Why?
Could $x_1$ be calculated in a different way?
Δ = np.float16(b**2-4*a*c)
x1 = -2*c/(np.sqrt(Δ) + b )
print ("x1 = ", x1, "Exact:", x1_exact)
print ("relative Error: {:.4f}%".format((x1-x1_exact)/x1*100))
x1 = -0.01333680646001563 Exact: -0.01333570454687738 relative Error: 0.0083%
Even for $\Delta$ with lower accuracy, the final result is almost 12000 times more accurate than before!
Sometimes it is enough to reformulate an equation to have better defined solution
But still: $$\frac{1}{3}x^2+\frac15x+\frac{3}{100}=0$$
a, b, c = 1/3, 0.2, 0.03
b**2-4*a*c # test Δ == 0 would fail
1.3877787807814457e-17
a, b, c = 100, 60, 9
b**2-4*a*c # Δ = 0 exactly
0
Other examples:
$$E_0=\frac{mc^2}{\sqrt{1-\frac{v^2}{c^2}}}-mc^2$$$$E_1 = \frac{mv^2}{\sqrt{1-\frac{v^2}{c^2}} + 1-\frac{v^2}{c^2}} $$$$E_2 = \frac12mv^2$$m, c, v, E = 1, 1e9, 0.5, np.empty(4) ;
E[0] = m*c**2/np.sqrt(1-v**2/c**2)-m*c**2
E[1] = m*v**2/(np.sqrt(1-v**2/c**2)+1-v**2/c**2)
E[2] = m*v**2/2
v=np.float128(v) # check what happens when v=0.1
E[3] = m*c**2/np.sqrt(1-v**2/c**2)-m*c**2
print("E = ", E)
E = [0. 0.125 0.125 0.1875]
T = np.float16
a, b, c = T(0.1), T(0.2), T(0.3)
print(1e16*(a+b-c))
a, b, c = T(0.1), T(0.2), T(a+b)
print(1e16*(a+b-c))
-2441406250000.0 0.0
"""Standard summation"""
data = np.array([1000]+1000*[0.1])
def naive_summation(data):
s = 0
for d in data:
s += d
return s
print(naive_summation(data))
1099.9999999999363
data2=np.copy(data) # data2 = data is just an alias to data (reference not a real copy)
data2.sort()
print(naive_summation(data2)) # better but not accurate
print(np.sum(data)) # standard numpy summation
1099.9999999999986 1100.0000000000002
"""Kahan Summation"""
def kahn_summation(data):
s, c = 0, 0
for d in data:
y = d-c
t = s+y
c = (t-s)-y
s = t
return s
print(kahn_summation(data))
1100.0
Question: how much a numerical procedure $f(x)$ multiples an input error $\Delta x$
Absolute condition number: how much an absolute error can be increased
Relative condition number:
$$\text{cond}(f)=\lim_{\varepsilon \rightarrow 0} \sup_{\|\Delta \| x \leq \varepsilon} \frac{\|\Delta f(x)\| / \|f(x)\|}{\|\Delta x\| / \|x\|}$$Condition number is a property of a given problem not the method.
Example from previous lecture $$\scriptsize I_n=\frac{1}{e}\int_0^1e^xx^n\;dx$$
generates the following conditions (decreasing but positive sequence)
$$\scriptsize 0<I_n<I_{n-1}$$and
$$\scriptsize I_0=1-\frac{1}{e}\qquad I_n=1-nI_{n-1}$$Suppose that we calculate $I_0$ with $\Delta I_0$ error
So the error of $I_n$ is
$$\Delta I_{n}=(-n)(-(n-1))(-(n-2))\ldots(-2)(-1)\,\Delta I_{0}=\Delta I_{n}=(-1)^nn!\,\Delta I_{0}$$So the condition number for calculation of $I_n$ is $n!$. Even if $\Delta I_0=10^{-16}$ we have $\Delta I_{20}\approx 240$ with $\text{cond}>2.4\cdot10^{18}$ since $I_n<I_0$
def F(x):
return x**4-2
import scipy.optimize as sc
x = sc.fsolve(F, 1)
print(x[0], 2**0.25)
1.189207115002721 1.189207115002721
Python (or more precisely scipy) provides a powerful solver
It is enough to provide a function definition and a starting point
But
Suppose that we have a continous function $f:\,\mathbb{R}\ni x\to\mathbb{R}$
Task: find a solution $f(x)=0$.
Two basic strategies:
Bisection method
Note: count $f$ only once for each point
make_plot()
In each step we reduce twice the legth of the interval in which the root resides
$\epsilon_n=\frac{1}{2}\epsilon_{n-1}=2^{-n}(b-a)$
Example of a linear convergence
$\epsilon_n\sim\epsilon^{\color{red}1}_{n-1}$
Why bisection and not decasection?
make_Newton_plot(0.6, f, df)
[0.6 1.60740741 1.26575079 1.07259388 1.0070429 ]
Newton method can converge very fast $x_n=x_*+\epsilon_n: f(x_*)=0$: $$\small (x_{n+1}-x_n)f'(x_n)=-f(x_n)$$ Expanding in Taylor series up to $\mathcal{O}(\epsilon^2)$:
$$\small (\epsilon_{n+1}-\epsilon_n)(f'(x_*)+\epsilon_n f''(x_*))=-f(x_*)-\epsilon_n f'(x_*)-\frac12\epsilon_n^2 f''(x_*)$$The convergence can be quadratic $$\small \epsilon_{n+1} = \frac{f''(x_*)}{2f'(x_*)}\epsilon_n^{\color{red}2}$$ provided that $f'(x_*)\neq 0$.
What happens when $f(x_*)=0$: exapmle: $f(x)=|x|^k$.
$$x_{n+1}=x_n-\frac{1}{k}x_n=\left(1-\frac1k\right)x_n=q\,x_n$$This is a geometric series which converges linearly $|q|<1$ or can be even divergent $|q|>1$.
Interesting case $q=-1$ for $k=1/2$
$$x_{n+1}=-x_n \qquad \text{two-cycle} $$make_Newton_plot(0.6, lambda x: np.sqrt(np.abs(x)), \
lambda x: 0.5*np.sign(x)/np.sqrt(np.abs(x)), \
left=-1, right=1, bottom =-0.5, top=1)
[ 0.6 -0.6 0.6 -0.6 0.6]
Newton method:
Similar methods:
Note: both Newton and Hayley developed their methods for celestial mechanical problems to solve Kepler's problem (finding E) $$M=E-e\sin E$$ where $M$ - mean anomaly, $E$ - eccentric anomaly, and $e$ - eccentricity.
Other high order methods?
Take parabola from three points $(x_n,y_n), (x_{n-1},y_{n-1}), (x_{n-2},y_{n-2})$ instead of a straight line?
Yes, but two solutions for $Q(x)=0$ (or one or none or infinitely many)
Muller's method - some more stable algorithm
Better: inverse interpolation
$$\scriptsize P(y)=\frac{(y-y_{n-1})(y-y_{n-2})x_n}{(y_n-y_{n-1})(y_n-y_{n-2})} + \frac{(y-y_{n})(y-y_{n-2})x_{n-1}}{(y_{n-1}-y_{n})(y_{n-1}-y_{n-2})} + \frac{(y-x_{n})(y-x_{n-1})x_{n-2}}{(y_{n-2}-y_{n})(y_{n-2}-y_{n-1})}$$$P(y)$ interpolates the inverse of $y=f(x)$ i.e. $x=f^{-1}(y)\approx P(y)$ so $x_{n+1}=P(y=0)$ gives another approximation to zero.
Polynomial equations $Q(x)=0$ can make use of some better techniques:
Use the scipy solver to solve and your own methods to solve the following equations
Find a solution of $z^3+1=0$ using Newton's method starting from complex snumbers.
Yield to Maturty.
YTM is a solution to an equation for $y$:
$$P=\frac{F}{(1+y)^T}+\sum_{i=1}^N\frac{C_i}{(1+y)^i}$$
Write a program to solve for $y$ for given bond price $P$, face value $F$ and coupons $C\in \mathbb{R}^T$ (narray), for e\xample choose the following values: $$P=\$89,\quad F=\$100, \quad C_i=\$3,\quad T=29$$