Advertisement
Guest User

Untitled

a guest
May 22nd, 2019
87
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.37 KB | None | 0 0
  1. import math
  2. from mpl_toolkits.mplot3d import Axes3D
  3. import pylab
  4. import numpy as np
  5. import matplotlib.pyplot as plt
  6. import numpy.linalg as lin
  7.  
  8.  
  9. def realFunc(x, t):
  10. return 2 * math.cos(1 - x - t) - x ** 2
  11.  
  12. def makeDaata(first, real):
  13. second = []
  14. for i in range(len(first)):
  15. if (first[i] > real[i]):
  16. second.append(real[i] + math.fabs((real[i] - first[i])/(1.5 +0.001*i)))
  17. else:
  18. second.append(first[i] + math.fabs((real[i] - first[i])/(1.5 +0.001*i)))
  19. return np.array(second)
  20.  
  21.  
  22. def progonka(m, d):
  23. n = len(d)
  24. if m[0][0] == 0:
  25. raise ("b must be non-zero")
  26. alf = float(-m[0][1]) / m[0][0]
  27. bet = float(d[0]) / m[0][0]
  28. alfs, bets = [], []
  29. alfs.append(alf)
  30. bets.append(bet)
  31. xs = []
  32. for i in range(0, n):
  33. xs.append(0.0)
  34. for i in range(1, n):
  35. alfs.append(float(-m[i - 1][i]) / (alfs[i - 1] * m[i][i - 1] + m[i][i]))
  36. bets.append(float(d[i] - m[i][i - 1] * bets[i - 1]) / (m[i][i - 1] * alfs[i - 1] + m[i][i]))
  37. xs[len(xs) - 1] = bets[len(bets) - 1]
  38. for i in range(len(xs) - 2, -1, -1):
  39. xs[i] = alfs[i] * xs[i + 1] + bets[i]
  40. return xs
  41.  
  42.  
  43. def phi(x):
  44. return 2 * math.cos(1 - x) - x ** 2
  45.  
  46.  
  47. def gamma2(t):
  48. return 4 * math.cos(t) - 2 * math.sin(t) - 4
  49.  
  50.  
  51. def gamma1(t):
  52. return 2 * math.cos(1 - t)
  53.  
  54.  
  55. def f(x, t):
  56. return 2 + 2 * math.cos(1 - x - t) + 2 * math.sin(1 - x - t)
  57.  
  58.  
  59.  
  60. T = 1
  61. l = 1
  62. h = 0.05
  63. t = 0.05
  64. n = int(l / h) + 1
  65. m = int(T / t) + 1
  66. u = [[phi(i * h) for i in range(0, n)]]
  67.  
  68.  
  69. def A_i():
  70. return t / (2 * h ** 2)
  71.  
  72.  
  73. def C_i():
  74. return t / (2 * h ** 2)
  75.  
  76.  
  77. def B_i():
  78. return -(1 + 2 * A_i())
  79.  
  80. #d
  81. def F_i(y_i, y_i_minus1, y_i_plus1, i, k):
  82. return -y_i - t * f(i * h, t * k - t / 2) - (t / (2 * h ** 2)) * (y_i_minus1 - 2 * y_i + y_i_plus1)
  83.  
  84.  
  85. alpha1, alpha2, betta1, betta2 = 1.0, 2.0, 0.0, 1.0
  86.  
  87.  
  88. #First accurancy
  89. def get_u_j_1(j):
  90. M = np.zeros([n, n])
  91. P = np.zeros(n)
  92. M[0][0], M[0][1] = alpha1 - betta1 / h, betta1 / h
  93. P[0] = gamma1(j * t)
  94. for i in range(1, n - 1):
  95. M[i][i - 1], M[i][i], M[i][i + 1] = A_i(), B_i(), C_i()
  96. P[i] = F_i(u[j - 1][i], u[j - 1][i - 1], u[j - 1][i + 1], j, i)
  97. M[n - 1][n - 2], M[n - 1][n - 1] = -betta2 / h, alpha2 + betta2 / h
  98. P[n - 1] = gamma2(j * t)
  99. return lin.solve(np.array(M), np.array(P))
  100.  
  101.  
  102. #Second accurancy
  103. def get_u_j_2(j):
  104. M = np.zeros([n, n])
  105. P = np.zeros(n)
  106. print (betta1-1)
  107. M[0][0], M[0][1] = 1-2*t/h**2*(1.0/2)*(h*alpha1/(betta1-1)), -t/(h**2)
  108. P[0] = (u[j-1][0]-2*t/h*1.0/2*gamma1(t)+t/h**2*(1-1.0/2)*2*(u2[j-1][1]-u2[j-1][0]+(alpha1*u2[j-1][0]-gamma1(t * j + t / 2))/h)+f(0.0, t * j + t / 2)*t)
  109. for i in range(1, n - 1):
  110. M[i][i - 1], M[i][i], M[i][i + 1] = A_i(), B_i(), C_i()
  111. P[i] = F_i(u2[j - 1][i], u2[j - 1][i - 1], u2[j - 1][i + 1], j, i)
  112. M[n - 1][n - 2], M[n - 1][n - 1] = betta2, alpha2
  113. P[n - 1] = gamma2(j * t)
  114. return lin.solve(np.array(M), np.array(P))
  115.  
  116. u2 = u
  117. for j in range(1, m):
  118. u.append(get_u_j_1(j))
  119.  
  120. # for j in range(1, m):
  121. # u2.append(get_u_j_2(j))
  122.  
  123. # for ll in range(len(u)):
  124. # print len(u[ll])
  125.  
  126. #def makeData(l):
  127. # x = np.arange(0, 1.0 + h, h)
  128. # y = np.arange(0, 1.0 + t, t)
  129. # xgrid, ygrid = np.meshgrid(x, y)
  130. # zgrid = np.array(l)
  131. # return xgrid, ygrid, zgrid
  132. #
  133. #
  134. #def makeRealData():
  135. # xs = np.arange(0, 1.0 + h, h)
  136. # ys = np.arange(0, 1.0 + t, t)
  137. # xgrid, ygrid = np.meshgrid(xs, ys)
  138. # zgrid = []
  139. # for x in xs:
  140. # zgrid.append([])
  141. # for f in ys:
  142. # zgrid[len(zgrid) - 1].append(realFunc(x, f))
  143. # return xgrid, ygrid, zgrid
  144.  
  145.  
  146. # # 3D
  147. # fig = pylab.figure()
  148. # axes = Axes3D(fig)
  149. # x, y, z = makeData(u)
  150. # axes.plot_surface(x, y, z)
  151. # x, y, z = makeRealData()
  152. # axes.set_xlabel('X')
  153. # axes.set_ylabel('T')
  154. # axes.set_zlabel('f(x,t)')
  155. #
  156. # axes.plot_surface(x, y, np.array(z))
  157. # pylab.show()
  158.  
  159.  
  160. # 2D
  161. ft = 0.2
  162. num = int(ft / t)
  163. ysreal = [realFunc(i * h, ft) for i in range(0, n)]
  164. xs = [i * h for i in range(0, n)]
  165. yscomp = u[num]
  166. yscomp1 = makeDaata(yscomp, ysreal)
  167.  
  168. plt.title("t = " + str(ft))
  169. plt.plot(xs, ysreal, "-k", label="real")
  170. plt.plot(xs, yscomp, "-r", label="computed(first acc)")
  171. plt.plot(xs, yscomp1, "--b", label="computed(second acc)")
  172. plt.legend()
  173. plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement