Advertisement
Guest User

Untitled

a guest
May 21st, 2019
94
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.37 KB | None | 0 0
  1. import math
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  4. import numpy
  5. from mpl_toolkits.mplot3d import Axes3D
  6. import pylab
  7.  
  8. def func1(j,h):
  9.     x = j*h
  10.     return x**2
  11.  
  12. def func2(j,h,t):
  13.     x = j*h
  14.     return x*t+x**2
  15.  
  16. def gamma1(k,t):
  17.     t = t*k
  18.     return -t
  19.  
  20. def gamma2(k,t):
  21.     t = k*t
  22.     return 2+t/((math.cosh(t)))**2
  23.  
  24. def realfunc(j,k,h,t):
  25.     x = j*h
  26.     t = k*t
  27.     return x**2+(math.tanh(t*x))
  28.  
  29. def  square(u,r):
  30.     x = 0.0
  31.     for  k in range (len(u)):
  32.         for j in range (len(u[k])):
  33.             x+=(u[k][j]-r[k][j])**2
  34.     return x
  35. def f(j,k,h,t):
  36.     x = h*j
  37.     t1 = t*j
  38.     f = -1 +((t1)**2-2*(x**2))*(math.tanh(x*t1)/(math.cosh(x*t1))**2)
  39.     return f
  40.  
  41.  
  42. l = 1
  43. T = 1
  44. h = 0.05
  45. t = 0.05
  46. n = int(l/h)
  47. m = int(T/t)
  48. g = (t**2)/(h**2)/2
  49. u =[
  50.     [],[]
  51. ]
  52.  
  53. for j in range(0,n):
  54.     u[0].append(func1(j,h))
  55.     u[1].append(func2(j,h,t))
  56.  
  57. for k in range(2,m):
  58.     u.append([])
  59.     for j in range(1,n-1):
  60.         u[k].append((g**2)*u[k-1][j-1]+(2 - 2*(g**2))*u[k-1][j]+(g**2)*u[k-1][j+1]-u[k-2][j])
  61.     u[k].insert(0,(1/(1+h)*u[k][0]+h*gamma1(k,t)/(h+1)))
  62.     u[k].append(u[k][n-2]+h*gamma2(k,t))
  63.  
  64. u2 = [
  65.     [],[]
  66. ]
  67.  
  68. for j in range(0,m):
  69.     u2[0].append(func1(j,h))
  70.  
  71. for j in range(1,m-1):
  72.     u2[1].append(u2[0][j]+t*func2(j,h,t)/2+f(j,0,h,t)*(t**2)/2+(u2[0][j+1]-2*u2[0][j]+u2[0][j-1])*(t**2)/(4*(h**2)))
  73.  
  74. u2[1].insert(0,((4*u2[1][0]-u2[1][1]-gamma1(0,t)*2*h)/3))
  75. u2[1].append(2*h*gamma2(m,t)/3+4*u[1][m-2]/3-(u[1][m-3])/3)
  76.  
  77. for k in range(2,m):
  78.     u2.append([])
  79.     for j in range(1,n-1):
  80.         u2[k].append((g**2)*u2[k-1][j-1]+(2 - 2*(g**2))*u2[k-1][j]+(g**2)*u2[k-1][j+1]-u2[k-2][j])
  81.     u2[k].insert(0,(4*u2[k][0]-u2[k][1]+2*gamma1(k,t)*h)/(3+2*h))
  82.     u2[k].append((2*h*gamma2(k,t)-u2[k][n-3]+4*u2[k][n-2])/3)
  83.  
  84. u0 = []
  85.  
  86. for k in range(0,n):
  87.     u0.append([])
  88.     for j in range(0,m):
  89.         u0[k].append(realfunc(j,k,h,t))
  90.  
  91.  
  92. print(np.matrix(u0))
  93. print('jopa')
  94. print(np.matrix(u))
  95. print('pisya')
  96. print(np.matrix(u2))
  97. print('asd',square(u,u0))
  98. print('asd',square(u2,u0))
  99.  
  100.  
  101.  
  102.  
  103. def makeData ():
  104.     x = numpy.arange(0,1,0.05)
  105.     y = numpy.arange(0,1,0.05)
  106.     xgrid, ygrid = numpy.meshgrid(x, y)
  107.     zgrid = numpy.array(u)
  108.     return xgrid, ygrid, zgrid
  109.  
  110. def makeDataReal():
  111.     x = numpy.arange(0,1,0.05)
  112.     y = numpy.arange(0,1,0.05)
  113.     xgrid, ygrid = numpy.meshgrid(x, y)
  114.     zgrid = numpy.array(u0)
  115.     return xgrid, ygrid, zgrid
  116.  
  117. def makeDataRealSection(t):
  118.     x = numpy.arange(0,1,0.05)
  119.     zgrid = numpy.array(u0[t])
  120.     return x, zgrid
  121.  
  122. def makeData2Section(t):
  123.     x = numpy.arange(0,1,0.05)
  124.     zgrid = numpy.array(u2[t])
  125.     return x, zgrid
  126.  
  127. def makeDataSection(t):
  128.     x = numpy.arange(0,1,0.05)
  129.     zgrid = numpy.array(u[t])
  130.     return x, zgrid
  131.  
  132.  
  133. def makeData2():
  134.     x = numpy.arange(0,1,0.05)
  135.     y = numpy.arange(0,1,0.05)
  136.     xgrid, ygrid = numpy.meshgrid(x, y)
  137.     zgrid = numpy.array(u2)
  138.     return xgrid, ygrid, zgrid
  139.  
  140. #1
  141. fig = pylab.figure()
  142. axes  = Axes3D(fig)
  143. x, y, z = makeData()
  144. axes.plot_surface(x,y,z)
  145. axes.set_title("№1 orange - real , blue - numerical")
  146. x, y, z = makeDataReal()
  147. axes.set_xlabel('X')
  148. axes.set_ylabel('T')
  149. axes.set_zlabel('Z')
  150. axes.plot_surface(x,y,z)
  151.  
  152. #2
  153. fig = pylab.figure()
  154. axes  = Axes3D(fig)
  155. x, y, z = makeData2()
  156. axes.plot_surface(x,y,z)
  157. axes.set_title("№2 orange - real , blue - numerical")
  158. x, y, z = makeDataReal()
  159. axes.set_xlabel('X')
  160. axes.set_ylabel('T')
  161. axes.set_zlabel('Z')
  162. axes.plot_surface(x,y,z)
  163. pylab.show()
  164.  
  165. #3 Section t = 0.5
  166. def generate2dData(t):
  167.     x1, y1 = makeDataSection(t)
  168.     x2,  y2 = makeData2Section(t)
  169.     x3, y3 = makeDataRealSection(t)
  170.     return x1,y1,x2,y2,x3,y3
  171.  
  172. def twoDemensionalGraphic(x1,y1,x2,y2,x3,y3,t):
  173.     plt.xlabel(r'$x$')
  174.     plt.ylabel(r'$x$')
  175.     plt.title("t =" + str(t))
  176.     plt.plot(x1, y1,"-k", label = "Первая точность")
  177.     plt.plot(x2, y2, "-r", label ="Вторая точность")
  178.     plt.plot(x3,y3, "-b", label="Реальное значение")
  179.     plt.legend()
  180.     plt.grid(True)
  181.     plt.show()
  182.  
  183. # меняй s(типо t), при котором будет строиться 2d график
  184. s = 0.1
  185. k = 0
  186. h = 0.0
  187. while True:
  188.     h += t
  189.     if (h >= s):
  190.         break
  191.     k += 1
  192. print(k,m)
  193. if (k>= m):
  194.     print("too big t")
  195. else:
  196.     x1,y1,x2,y2,x3,y3 = generate2dData(k)
  197.     twoDemensionalGraphic(x1,y1,x2,y2,x3,y3,s)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement