Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Project: WISB251
Views: 81
Kernel: Python 3 (Ubuntu Linux)
import numpy as np import matplotlib.pyplot as plt

Hoofdstuk 1 & 2

Benaderen van de afgeleide - nauwkeurigheid

x = 1.2; derivative = np.cos(x) for h in np.flip(np.logspace(-20,0,21),0): approximation1 = (np.sin(x + h) - np.sin(x))/h approximation2 = (np.sin(x + h) - np.sin(x - h))/(2*h) error1 = abs(derivative - approximation1) error2 = abs(derivative - approximation2) print(h, error1, error2)
1.0 0.485900436624 0.0574442179644 0.1 0.047166759977 0.000603627697885 0.01 0.00466619586072 6.0392657113e-06 0.001 0.000466079897112 6.03930143672e-08 0.0001 4.66025575812e-05 6.03618710571e-10 1e-05 4.66019133444e-06 2.43299824731e-12 1e-06 4.65968587493e-07 7.98405785929e-12 1e-07 4.61932618823e-08 1.19006360322e-10 1e-08 4.36105151991e-10 4.36105151991e-10 1e-09 5.59472563832e-08 4.36105151991e-10 1e-10 1.66969558846e-07 3.88141953467e-07 1e-11 7.93853073122e-06 2.3874156081e-06 1e-12 0.00013006306344 7.45519122087e-05 1e-13 0.000685174575753 0.00013006306344 1e-14 0.00401584364963 0.0015352714735 1e-15 0.0817314553734 0.0262203041421 1e-16 0.362357754477 0.362357754477 1e-17 0.362357754477 0.362357754477 1e-18 0.362357754477 0.362357754477 1e-19 0.362357754477 0.362357754477 1e-20 0.362357754477 0.362357754477

Evalueren van polynomen - efficientie, nauwkeurigheid en robuustheid

f1 = lambda x: x**7 - 7*x**6 + 21*x**5 - 35*x**4 + 35*x**3 - 21*x**2 + 7*x - 1 f2 = lambda x: -1 + x*(7 + x*(-21 + x*(35 + x*(-35 + x*(21 + x*(-7 + x)))))) f3 = lambda x: (x-1)**7 x = np.linspace(0.99,1.01,100) plt.plot(x,f1(x)) plt.plot(x,f2(x)) plt.plot(x,f3(x)) plt.show()
Image in a Jupyter notebook

Benaderen van 2\sqrt{2} - efficientie, nauwkeurigheid en robuustheid

xk = 100 n = 10 for k in range(n): xk = (xk + 2/xk)/2 ek = abs(xk - np.sqrt(2)) print(k, xk, ek)
0 50.01 48.5957864376 1 25.024996000799838 23.6107824384 2 12.552458046745903 11.1382444844 3 6.35589469493114 4.94168113256 4 3.335281609280434 1.92106804691 5 1.967465562231149 0.553251999858 6 1.4920008896897232 0.0777873273166 7 1.4162413320389438 0.00202776966585 8 1.4142150140500531 1.45167695798e-06 9 1.41421356237384 7.44959649523e-13

Oplossen van een stelsel vergelijkingen - nauwkeurigheid, robuustheid

import Gauss_elimination as GE A = np.array([[0.002, 1.231, 2.471],[1.196, 3.165, 2.543],[1.475, 4.271, 2.142]],dtype='float32') b = np.array([3.704, 6.904, 7.888],dtype='float32') xt = np.array([1,1,1],dtype='float32') GE.forward(A,b) x = GE.backward(A,b) print(np.linalg.norm(A.dot(x) - b)/np.linalg.norm(b)) print(np.linalg.norm(xt - x)/np.linalg.norm(xt))
5.39875e-11 0.000250739
n = 100 A = np.eye(n) - np.tri(n,k=-1) A[:,n-1]=1 xt = np.random.randn(n) b = A.dot(xt) GE.forward(A,b) x = GE.backward(A,b) print(np.linalg.norm(A.dot(x) - b)/np.linalg.norm(b)) print(np.linalg.norm(xt - x)/np.linalg.norm(xt))
0.0 0.708478719517

Benaderen van een integraal - ophopen van afrondfouten

n = 30 y = np.zeros(n) y[0] = np.log(11) - np.log(10) for k in range(1,n): y[k] = 1.0/n - 10*y[k-1] print(y)
[ 9.53101798e-02 -9.19768465e-01 9.23101798e+00 -9.22768465e+01 9.22801798e+02 -9.22798465e+03 9.22798798e+04 -9.22798765e+05 9.22798768e+06 -9.22798768e+07 9.22798768e+08 -9.22798768e+09 9.22798768e+10 -9.22798768e+11 9.22798768e+12 -9.22798768e+13 9.22798768e+14 -9.22798768e+15 9.22798768e+16 -9.22798768e+17 9.22798768e+18 -9.22798768e+19 9.22798768e+20 -9.22798768e+21 9.22798768e+22 -9.22798768e+23 9.22798768e+24 -9.22798768e+25 9.22798768e+26 -9.22798768e+27]

Oplossen van een differentiaalvergelijking - nauwkeurigheid en stabiliteit

a = -10 T = 10 dt = 1 n = int(T/dt) t = np.linspace(0,T,n) y = np.zeros(n) y[0] = 1 for k in range(n-1): y[k+1] = y[k] + dt*a*y[k] plt.plot(t,y,t,np.exp(a*t))
[<matplotlib.lines.Line2D at 0x7f8538c63f98>, <matplotlib.lines.Line2D at 0x7f8538c6a9e8>]
Image in a Jupyter notebook

Opdelen van de reeele as

x = [] beta = 2 t = 3 L = -2 U = 3 for e in range(L,U+1): for d in range(1,beta**(t)): x.append(d*beta**(e-(t-1))) x.sort() print(x) plt.plot(x,np.zeros(len(x)),'.') plt.xlim((0,15))
[0.0625, 0.125, 0.125, 0.1875, 0.25, 0.25, 0.25, 0.3125, 0.375, 0.375, 0.4375, 0.5, 0.5, 0.5, 0.625, 0.75, 0.75, 0.875, 1.0, 1.0, 1, 1.25, 1.5, 1.5, 1.75, 2.0, 2, 2, 2.5, 3.0, 3, 3.5, 4, 4, 5, 6, 6, 7, 8, 10, 12, 14]
(0, 15)
Image in a Jupyter notebook

Rekenen met floating point getallen

We kunnen met Python met minder nauwkeurigheid rekenen door float16 of float32 te gebruiken

eta = 2**-10 print('eta = ', eta) x = np.random.randn() print(abs(x - np.float16(x))/abs(x) < eta)
eta = 0.0009765625 True

Wat kan er mis gaan?

# Verlies van significantie x = np.float16(12541) y = np.float16(2.1421e-4) print(float(x)+float(y), x+y) print(y,(x+y)-x)
12544.000214219093 12544.0 0.00021422 0.0
# gebeurt niet alleen met 16-bits float! x = 12541 y = 2.1421e-4 print(y,(x+y)-x,(x-x)+y)
0.00021421 0.00021421000019472558 0.00021421
# Delen door kleine getallen x = np.random.randn() y = np.random.randn()*1e-4 print(abs((x / y))) print(abs((x / y) - np.float16(np.float16(x) / np.float16(y)))) print(abs((x / y) - np.float16(np.float16(x) / np.float16(y))) / abs((x / y)))
21491.926601766816 19.9266017668 0.000927166844371
# Cancellation x = 1.43145542 y = x + 1e-4 print(abs((x - y) - np.float16(np.float16(x) - np.float16(y))) / abs((x - y)))
1.0
# Overflow x = 1.43e4 y = 2.32 print(np.sqrt(np.float16(x**2) + np.float16(y**2))) print(np.float16(np.sqrt(np.float16(1) + np.float16((y/x)**2))*np.float16(x)))
inf 14304.0
# Underflow x = 1.321e-4 y = 2.43e4 print(np.float16(x/y))
0.0

De volgorde is ook belangrijk!

1/3
0.3333333333333333
n = 10000 result1 = np.float16(0) for k in range(1,n): result1 = result1 + np.float16(1/k) result2 = np.float16(0) for k in range(n+1,0,-1): result2 = result2 + np.float16(1/k) result3 = np.float128(0.0) for k in range(1,n): result3 = result2 + 1/k print(result1, result2, result3)
7.0859 9.7969 9.79697501