Numerical Methods¶
Lecture 3a: Introduction to numerical methods¶
by: Tomasz Romańczukiewicz
web: th.if.uj.edu.pl/~trom/
rm B-2-03
Outline¶
- Solutions
Implementaion of the Newton's method
$$x\leftarrow x-\frac{f(x)}{f'(x)}$$for
$$f(x)=x-\cos(x),\qquad f'(x)=1+\sin(x)$$In [1]:
import numpy as np
x = 0
for n in range(10):
f, df = x-np.cos(x), 1+np.sin(x)
x -= f/df
print(x)
We can define a function returning both the value and the derivative
In [2]:
def function(x):
return x-np.cos(x), 1+np.sin(x)
def newton_method(F, x, N):
for n in range(N):
f, df = F(x)
x -= f/df
print(x)
newton_method(function, 1, 10)
In [3]:
from scipy.optimize import fsolve
help(fsolve)
But functions can have more arguments.
fsolve(func, x0, args=() ...)
allows arguments as a tuple.
In [5]:
def test_function(func, argc=()):
func(*argc) # execute func with expanded tuple argc
def g(a, b, c):
print("g: a = ", a, ", b = ", b, ", c =", c)
my_args=(666, "Hello world", 3.1416) # arguments passed as a tuple
test_function(g, my_args)
In [6]:
def function(x, a, b):
return a+b*x-np.cos(x), b+np.sin(x)
def newton_method(F, x0, N, args=()):
x = x0
for n in range(N):
f, df = F(x, *args)
x -= f/df
print(x)
newton_method(function, 1, 10, args=(0,1))
In [7]:
""" Add more default arguments """
def newton_method(F, x, N=100, args=(), xtol=1e-8):
for n in range(N):
f, df = F(x, *args)
dx = f/df
x -= dx
if np.abs(dx)<xtol: break
return (x, f, dx, n)
result = newton_method(function, 1, 10, args=(0,1))
print(result)
In [65]:
"""Returning a dictionary for better understanding"""
def newton_method(F, x, N=100, args=(), xtol=1e-8):
for n in range(N):
f, df = F(x, *args)
dx = f/df
x -= dx
if np.abs(dx)<xtol: break
return {'x': x, 'lastf': F(x,*args)[0], 'last_step': dx, 'iterations': n}
result = newton_method(function, 1, 10, args=(0,1))
print(result)
In [66]:
result = newton_method(lambda x: (x**4-2, 4*x**3), 1)
print(result)
In [69]:
"""Using sympy to calculate derivatives"""
from sympy.abc import x
import sympy as sp
f = x**4-2
lam_f = sp.lambdify(x, (f, sp.diff(f)))
print("", newton_method(lam_f, 1))
In [10]:
%matplotlib inline
import matplotlib.pyplot as plt
N = 1000
x = np.linspace(-1,1, N)
y = np.cos(10*np.arccos(x))
plt.plot(x,y)
plt.show()
In [11]:
z = np.empty([N, 8])
for n in range(8):
plt.plot(x, np.cos(n*np.arccos(x)))
plt.show()
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$$
In [12]:
def P_python(y, F, T, C):
z = 1+y
P = F/z**T
for i in range(1, T+1):
P += C[i-1]/z**i
return P-89
F, T, C = 100, 29, [3]*29
print( P_python(0.02, F, T, C) )
In [13]:
import matplotlib.pyplot as plt
import numpy as np
Y = np.linspace(0, .1, 100)
plt.plot(Y, P_python(Y, F, T, C))
plt.grid(True)
plt.show()
In [14]:
def P_Cstyle(y, F, T, C):
z = 1+y
zi = 1
P = 0
for i in range(1, T+1):
zi*=z
P += C[i-1]/zi
return P+F/zi-89
print( P_Cstyle(0.02, F, T, C) )
In [15]:
"""Perhaps better with vectorization"""
def P_broadcasting(y, F, T, C):
z = 1+y
i = np.arange(1,T+1)
P = np.sum(C/z**i) + F/z**T
return P-89
F, T, C = 100, 29, np.array([3]*29)
print( P_broadcasting(0.02, F, T, C) )
In [50]:
import numba as nb
@nb.jit(nopython=True, cache=True)
def P_numba(y, F, T, C):
z = 1+y
i = np.arange(1,T+1)
P = np.sum(C/z**i) + F/z**T
return P-89
F, T, C = 100, 29, np.array([3]*29)
print( P_broadcasting(0.02, F, T, C) )
In [49]:
@nb.jit(nopython=True, cache=True)
def P_numba2(y, F, T, C):
z = 1+y
zi = 1
P = 0
for i in range(1, T+1):
zi*=z
P += C[i-1]/zi
return P+F/zi-89
F, T, C = 100, 29, np.array([3]*29)
print( P_numba2(0.02, F, T, C) )
In [47]:
import time
def test_function(function, ref_time=0):
T=20000
F, C = 0, np.array([0.1]*T)
t1 = time.time()
function(0.02, F, T, C)
t2 = time.time()
if ref_time==0: ref_time = (t2-t1)*1e3
return [ function.__name__, (t2-t1)*1e3, ref_time/(t2-t1)*1e-3, ""]
In [51]:
times = [test_function(P_python)]
for f in [P_Cstyle, P_broadcasting, P_numba, P_numba, P_numba2, P_numba2]:
times.append(test_function(f, times[0][1]))
times[-1][-1] = "Speed King"
import pandas as pd
pd.DataFrame(times, columns=["Function", "Time", "x faster then Python", ""])
Out[51]:
Question:
How long would it take to calculate $P$ using PySym rational numbers?
How long would it take to calculate $P$ using PySym rational numbers?
In [22]:
import numpy as np
import matplotlib.pyplot as plt
# %% Newton
def Newton_fractal_vectorize ( h,w, maxit=50 ):
"""Returns an image of the Newton fractal of size (h,w)."""
y,x = np.ogrid[ -2:2:h*1j, -2:2:w*1j ]
z = x+y*1j
for i in range(maxit):
#z -= (z**3-1)/(3*z**2)
z -= (z**3-2*z+2)/(3*z**2-2)
return np.angle(z+1j-1)
In [23]:
plt.figure(figsize=(10,10))
plt.imshow(Newton_fractal_vectorize(300,300))
plt.show()
