Chunk / Numerical

Download as
Browse Code »

Commit 7d0297f3f8739a34b4b89b3ecc5deb3862e92d1a

Authored by Chunk
1 parent e5dbe416
Exists in master

6 finished.

Showing 8 changed files with 285 additions and 4 deletions   Show diff stats
chap3/__init__.pyc 0 → 100644
View file @7d0297f
No preview for this file type
chap3/cholesky.py
  Diff comments   View file @7d0297f
... ... @@ -42,6 +42,14 @@ def gen_b(n, x=None):
42 42 return np.dot(gen_hilbert(n), x)
43 43  
44 44  
  45 +def norm(l):
  46 + s = 0
  47 + for item in l:
  48 + if np.absolute(item) > s:
  49 + s = np.absolute(item)
  50 + return s
  51 +
  52 +
45 53 def cholesky(A):
46 54 n = len(A)
47 55 for j in range(n):
... ... @@ -57,7 +65,7 @@ def cholesky(A):
57 65 return np.array(A)
58 66  
59 67  
60   -def calculate(n):
  68 +def calculate0(n):
61 69 H = gen_hilbert(n)
62 70 b = gen_b(n)
63 71 L = cholesky(H)
... ... @@ -68,6 +76,19 @@ def calculate(n):
68 76 return x, b - np.dot(H, x), [x[i] - 1.0 for i in range(n)]
69 77  
70 78  
  79 +def calculate1(n):
  80 + H = gen_hilbert(n)
  81 + b = gen_b(n)
  82 + b_t = np.copy(b)
  83 + b_t[-1] += math.pow(10, -7)
  84 + L = cholesky(H)
  85 + y = calc_lower(L, b_t)
  86 + x = calc_upper(L.T, y)
  87 + # validate:
  88 + # print np.dot(H, x) - b
  89 + return x, b - np.dot(H, x), [x[i] - 1.0 for i in range(n)]
  90 +
  91 +
71 92 def test():
72 93 print gen_hilbert(3)
73 94 print gen_b(5)
... ... @@ -83,10 +104,30 @@ def test():
83 104 A = [[5, -1, -1], [-1, 3, -1], [-1, -1, 5]]
84 105 print cholesky(A)
85 106  
86   - x, r, delta = calculate(10)
87   - print "%s\n%s\n%s\n" % (x, r, delta)
  107 +
  108 +def test1():
  109 + print "[origin(n=10)]"
  110 + x, r, delta = calculate0(10)
  111 + print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
  112 + print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
  113 +
  114 + print "[disturbed(n=10)]"
  115 + x, r, delta = calculate1(10)
  116 + print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
  117 + print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
  118 +
  119 + print "[n=8]"
  120 + x, r, delta = calculate0(8)
  121 + print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
  122 + print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
  123 +
  124 + print "[n=12]"
  125 + x, r, delta = calculate0(12)
  126 + print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
  127 + print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
88 128  
89 129  
90 130 if __name__ == '__main__':
91   - test()
  131 + # test()
  132 + test1()
92 133 pass
... ...
chap3/cholesky.pyc 0 → 100644
View file @7d0297f
No preview for this file type
chap4/iterative.py 0 → 100644
  Diff comments   View file @7d0297f
... ... @@ -0,0 +1,87 @@
  1 +__author__ = 'chunk'
  2 +
  3 +import numpy as np
  4 +import math
  5 +
  6 +
  7 +def gen_hilbert(n):
  8 + return np.array([1.0 / (i + j + 1) for i in range(n) for j in range(n)]).reshape((n, -1))
  9 +
  10 +
  11 +def gen_b(n):
  12 + return np.array([1.0 / i for i in range(1, n + 1)])
  13 +
  14 +
  15 +def norm(l):
  16 + s = 0
  17 + for item in l:
  18 + if np.absolute(item) > s:
  19 + s = np.absolute(item)
  20 + return s
  21 +
  22 +
  23 +def jacobi(A, b, epsilon=0.0001, x=None):
  24 + n = len(A)
  25 + assert n == len(b)
  26 +
  27 + if x is None:
  28 + x = np.zeros(n, dtype=np.float32)
  29 + x_old = None
  30 +
  31 + D = np.diag(A)
  32 + R = A - np.diagflat(D)
  33 +
  34 + # for i in range(N):
  35 + iter = 0
  36 + while x_old is None or norm(x - x_old) >= epsilon:
  37 + iter += 1
  38 + x_old = np.copy(x)
  39 + x = (b - np.dot(R, x_old)) / D
  40 + print iter, x
  41 + return x
  42 +
  43 +
  44 +def SOR(A, b, omega, epsilon=0.0001, x=None):
  45 + n = len(A)
  46 + assert n == len(b)
  47 +
  48 + if x is None:
  49 + x = np.zeros(n, dtype=np.float32)
  50 + x_old = None
  51 + iter = 0
  52 + while x_old is None or norm(x - x_old) >= epsilon:
  53 + iter += 1
  54 + x_old = np.copy(x)
  55 + for i in range(n):
  56 + s = 0.0
  57 + for j in range(n):
  58 + if j != i:
  59 + s += A[i][j] * x[j]
  60 + x[i] = (1 - omega) * x[i] + omega * (b[i] - s) / A[i][i]
  61 +
  62 + # print iter, x
  63 + return iter, x
  64 +
  65 +
  66 +def test():
  67 + # A = [[10, 3, 1], [2, -10, 3], [1, 3, 10]]
  68 + # b = [14, -5, 14]
  69 + # print SOR(A, b, 1.25)
  70 + # print jacobi(A, b)
  71 +
  72 + n = 10
  73 + A = gen_hilbert(n)
  74 + b = gen_b(n)
  75 + # print jacobi(A, b)
  76 + print SOR(A, b, 1)
  77 +
  78 + for w in np.linspace(0.5, 1.5, num=11):
  79 + iter, x = SOR(A, b, w)
  80 + # print "w:%f\titer:%d\tx:%s\t" % (w, iter, x)
  81 + # print "norm(delta):%s\n" % (norm(np.array([1.0] + [0.0] * (n - 1)) - x))
  82 + print "w:%f\titer:%d\tnorm(delta):%s\n" % (
  83 + w, iter, norm(np.array([1.0] + [0.0] * (n - 1)) - x))
  84 +
  85 +
  86 +if __name__ == '__main__':
  87 + test()
... ...
chap5/__init__.py 0 → 100644
  Diff comments   View file @7d0297f
... ... @@ -0,0 +1 @@
  1 +__author__ = 'chunk'
... ...
chap5/power.py 0 → 100644
  Diff comments   View file @7d0297f
... ... @@ -0,0 +1,62 @@
  1 +__author__ = 'chunk'
  2 +
  3 +import numpy as np
  4 +import math
  5 +
  6 +
  7 +def norm(l):
  8 + s = 0
  9 + for item in l:
  10 + if np.absolute(item) > s:
  11 + s = np.absolute(item)
  12 + return s
  13 +
  14 +
  15 +def norm_vec(l):
  16 + s = 0
  17 + loc = 0
  18 + i = -1
  19 + for item in l:
  20 + i += 1
  21 + if np.absolute(item) > s:
  22 + s = np.absolute(item)
  23 + loc = i
  24 + return l[loc], np.array(l, dtype=np.float32) / l[loc]
  25 +
  26 +
  27 +def power_method(A, epsilon=0.00001, v=None):
  28 + n = len(A)
  29 + if v == None:
  30 + v = np.random.rand(n)
  31 + u = np.copy(v)
  32 + lmbda = None
  33 + lmbda_old = None
  34 +
  35 + iter = 0
  36 + while lmbda_old is None or np.absolute(lmbda - lmbda_old) >= epsilon:
  37 + iter += 1
  38 + lmbda_old = lmbda
  39 + v = np.dot(A, u)
  40 + lmbda, u = norm_vec(v)
  41 + print iter, lmbda
  42 +
  43 + return lmbda, u
  44 +
  45 +
  46 +def test():
  47 + ll = [0, 2, 0.5, 1, -3, 0.2]
  48 + ll2 = [i * 2 for i in ll]
  49 + print norm_vec(ll)
  50 + print norm_vec(ll2)
  51 +
  52 + A = [[3, 1], [1, 3]]
  53 + print power_method(A)
  54 +
  55 + A = [[5, -4, 1], [-4, 6, -4], [1, -4, 7]]
  56 + B = [[25, -41, 10, -6], [-41, 68, -17, 10], [10, -17, 5, -3], [-6, 10, -3, 2]]
  57 + print power_method(A)
  58 + print power_method(B)
  59 +
  60 +
  61 +if __name__ == '__main__':
  62 + test()
... ...
chap6/__init__.py 0 → 100644
  Diff comments   View file @7d0297f
... ... @@ -0,0 +1 @@
  1 +__author__ = 'chunk'
... ...
chap6/fitting.py 0 → 100644
  Diff comments   View file @7d0297f
... ... @@ -0,0 +1,89 @@
  1 +__author__ = 'chunk'
  2 +
  3 +import numpy as np
  4 +import math
  5 +import matplotlib.pyplot as plt
  6 +import seaborn as sns
  7 +
  8 +from chap3.cholesky import calc_upper, calc_lower, cholesky
  9 +
  10 +plt.ticklabel_format(style='sci', axis='both')
  11 +
  12 +
  13 +def fiiting(t, f, phi):
  14 + m, n = len(t), len(phi)
  15 + assert m == len(f)
  16 +
  17 + A = np.array([phi[j](t[i]) for i in range(m) for j in range(n)]).reshape(m, -1)
  18 + G = np.dot(A.T, A)
  19 + b = np.dot(A.T, f)
  20 + L = cholesky(G)
  21 + y = calc_lower(L, b)
  22 + x = calc_upper(L.T, y)
  23 + return x
  24 +
  25 +
  26 +def test0():
  27 + t0 = [1, 2, 3, 4, 5]
  28 + f0 = [4, 4.5, 6, 8, 8.5]
  29 + phi = [lambda x: 1, lambda x: x]
  30 + coef = fiiting(t0, f0, phi)
  31 + print coef
  32 +
  33 + fit = lambda x: coef[0] + coef[1] * x * x
  34 + t1 = np.linspace(1, 8, num=1000).tolist()
  35 + f1 = [fit(i) for i in t1]
  36 +
  37 + plt.scatter(t0, f0)
  38 + plt.plot(t1, f1)
  39 + plt.show()
  40 +
  41 +
  42 +def test():
  43 + # print np.linspace(1, 8, num=15).tolist()
  44 +
  45 +
  46 +
  47 + t = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0]
  48 + f = [33.40, 79.50, 122.65, 159.05, 189.15, 214.15, 238.65, 252.2, 267.55, 280.50, 296.65,
  49 + 301.65, 310.40,
  50 + 318.15, 325.15]
  51 + phi = [lambda x: 1, lambda x: x, lambda x: x * x]
  52 + coef = fiiting(t, f, phi)
  53 + print coef
  54 + fit = lambda x: coef[0] + coef[1] * x + coef[2] * x * x
  55 + t1 = np.linspace(1, 8, num=1000).tolist()
  56 + f1 = [fit(i) for i in t1]
  57 +
  58 +
  59 + ff = [np.log(i) for i in f]
  60 + print ff
  61 + phi2 = [lambda x: 1, lambda x: x]
  62 + coef2 = fiiting(t, ff, phi2)
  63 + print coef2
  64 + fit2 = lambda x: np.exp(coef2[0]) * np.exp(coef2[1] * x)
  65 + t2 = np.linspace(1, 8, num=1000).tolist()
  66 + f2 = [fit2(i) for i in t2]
  67 +
  68 + tt = [1.0/i for i in t]
  69 + ff = [np.log(i) for i in f]
  70 + print ff
  71 + phi2 = [lambda x: 1, lambda x: x]
  72 + coef2 = fiiting(tt, ff, phi2)
  73 + print coef2
  74 + fit2 = lambda x: np.exp(coef2[0]) * np.exp(coef2[1] * 1.0 / x)
  75 + t3 = np.linspace(1, 8, num=1000).tolist()
  76 + f3 = [fit2(i) for i in t3]
  77 +
  78 + plt.scatter(t, f)
  79 + plt.plot(t1, f1, label='polynomial')
  80 + plt.plot(t2, f2, label='exponential', linestyle='--')
  81 + plt.plot(t3, f3, label='exponential 1/t', linestyle='-')
  82 + plt.xlabel("t")
  83 + plt.ylabel("y")
  84 + plt.legend(loc=2)
  85 + plt.show()
  86 +
  87 +
  88 +if __name__ == '__main__':
  89 + test()
... ...