SHARE
TWEET

Untitled

a guest Oct 20th, 2019 82 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import math
  2.  
  3. N = 1000
  4.  
  5. EPS = 1e-4
  6.  
  7. ALPHA = 0.01
  8. BETA = 10
  9.  
  10.  
  11. def mul(bottom, middle, top, x):
  12.     ret = [0] * N
  13.  
  14.     for i in range(N):
  15.         ret[i] = middle[i] * x[i]
  16.         if i > 0:
  17.             ret[i] += bottom[i] * x[i - 1]
  18.         if i < N - 1:
  19.             ret[i] += top[i] * x[i + 1]
  20.  
  21.     return ret
  22.  
  23.  
  24. def diff(u, v):
  25.     return [u[i] - v[i] for i in range(N)]
  26.  
  27.  
  28. def norm(u):
  29.     return math.sqrt(sum(x * x for x in u))
  30.  
  31.  
  32. def universal_solver(bottom, middle, top, f, x, tau=None):
  33.     ret_e = []
  34.     inv = [1.0 / w for w in middle]
  35.     Ax = mul(bottom, middle, top, x)
  36.     r0 = norm(diff(Ax, f))
  37.     while True:
  38.         if tau is not None:
  39.             Ax = mul(bottom, middle, top, x)
  40.             Ax = [tau * w for w in Ax]
  41.         else:
  42.             Ax = mul([0] * N, inv, [0] * N, mul(bottom, middle, top, x))
  43.  
  44.         r = norm(diff(Ax, f))
  45.         e = r / r0
  46.  
  47.         ret_e.append(e)
  48.         if e < EPS:
  49.             break
  50.  
  51.         x = [x[i] + (f[i] - Ax[i]) for i in range(N)]
  52.  
  53.     return ret_e, x
  54.  
  55.  
  56. def main():
  57.     tau = 0.1651512666112257
  58.     bottom = [-1] * N
  59.     top = [-1] * N
  60.     middle = [2 + ALPHA] * N
  61.     middle[0] = 2 + BETA
  62.     x = [0] * N
  63.  
  64.     f = [1 if 494 <= i < 505 else 0 for i in range(N)]
  65.     f_tau = [tau * w for w in f]
  66.     inv_f = [1.0 / middle[i] if 494 <= i < 505 else 0 for i in range(N)]
  67.  
  68.     mti_e, mti_x = universal_solver(bottom, middle, top, f_tau, [0] * N, tau)
  69.     jacobi_e, jacobi_x = universal_solver(bottom, middle, top, inv_f, x)
  70.  
  71.     index_nev = [i + 1 for i in range(len(mti_e))]
  72.     plt.plot(index_nev, mti_e)
  73.     plt.show()
  74.  
  75.     index_solve = [i + 1 for i in range(len(mti_x))]
  76.     plt.plot(index_solve, mti_x)
  77.     plt.show()
  78.  
  79.  
  80. if __name__ == '__main__':
  81.     main()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top