• API
• FAQ
• Tools
• Archive
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 =  * 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( * N, inv,  * 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 = 2 + BETA
62.     x =  * 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,  * 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.

Top