Ledger Nano X - The secure hardware wallet
SHARE
TWEET

Untitled

a guest Apr 9th, 2020 146 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. from numpy import cos
  2. from numpy import sin
  3. from numpy import linspace
  4. import numpy as np
  5. import matplotlib.pyplot as plt
  6. import time
  7. from datetime import datetime
  8. from numpy import pi
  9. import copy
  10. start_time = datetime.now()
  11. from mpl_toolkits.mplot3d import Axes3D
  12. N10 = 10
  13. N20 = 20
  14. def exactsolution(x,y):
  15.     return cos(x**2+y)*sin(2*y+x)*pow(x**2+y**2+2*y+1,0.5)+sin(x**2+y)*cos(4*x+y)*pow(x+y+2,0.5)
  16. def border(number,x,y):#1 - фи1, 2 - фи2, 3- фи3, 4 -фи4
  17.     if (number ==1 ):
  18.         return exactsolution(0,y)
  19.     if (number == 2):
  20.         return exactsolution(1,y)
  21.     if (number == 3):
  22.         return exactsolution(x, 0)
  23.     if (number == 4):
  24.         return exactsolution(x, 1)
  25. def function(X,Y):#надо взять первую размерность массива
  26.     h1 = pow(10,-6)
  27.     h2 = pow(10,-6)/2
  28.     du_dx = ((exactsolution(X-h1,Y)-2*exactsolution(X,Y)+exactsolution(X+h1,Y))/(pow(h1,2)))
  29.     du_dy = ((exactsolution(X,Y-h2)-2*exactsolution(X,Y)+exactsolution(X,Y+h2))/(pow(h2,2)))
  30.     f = du_dx+du_dy
  31.     return f
  32. def progonka(A,B,C,D):#для уравнения A*u[j-1] + C*u[j] +B*u[j+1] = D[j]
  33.     N=len(D)
  34.     U = np.zeros(N)
  35.     P=np.zeros(N)
  36.     Q = np.zeros(N)
  37.     P[0] = -B[0]/C[0]
  38.     Q[0] = D[0]/C[0]
  39.     for j in range(1,N-1):
  40.         P[j] = -B[j]/(C[j]+A[j-1]*P[j-1])
  41.         Q[j] = (-A[j-1]*Q[j-1]+D[j])/(C[j]+A[j-1]*P[j-1])
  42.     last = N-1
  43.     U[last] = (-A[last-1]*Q[last-1]+D[last])/(C[last]+A[last-1]*P[last-1])
  44.     k=last
  45.     while(k!=0):
  46.         U[k-1]=P[k-1]*U[k]+Q[k-1]
  47.         k=k-1
  48.     return U
  49.  
  50. def solution(u,x,y,NX,NY):
  51.     l = len(x)
  52.     h1 = 1/NX
  53.     h2 = 1/NY
  54.     u_quad = fill_array(l,x,y)
  55.     u_half = u_quad.copy()#fill_array(l,x,y)
  56.     u_n1 =u_quad.copy()#fill_array(l,x,y)
  57.     Function = function(x,y)
  58.     right_vector1 =u_quad.copy()#fill_array(l,x,y)
  59.     right_vector2 =u_quad.copy()# fill_array(l,x,y)
  60.     for j in range(1,l-1):
  61.         for i in range(1,l-1):
  62.             right_vector1[i][j] = Function.copy()[i][j] - (2*u[i][j].copy() / tau)
  63.         #print(datetime.now() - start_time)
  64.         A = np.full(l-1,1/pow(h1,2),dtype=np.float64)
  65.         B = np.full(l-1,1/pow(h1,2),dtype=np.float64)
  66.         C = np.full(l,-(2/tau)-(2/pow(h1,2)),dtype=np.float64)
  67.         C[0] = 1.
  68.         B[0] = 0.
  69.         A[-1] = 0
  70.         C[-1] = 1.
  71.         u_quad[:,j] = progonka(A,B,C,right_vector1.copy()[:,j])
  72.         #print(datetime.now() - start_time)
  73.     for i in range(1,l-1):
  74.         for j in range(1,l-1):
  75.             right_vector2[i][j] = -(2*u_quad.copy()[i][j]/tau)
  76.         A1 = np.full(l-1, 1 / pow(h2, 2))
  77.         B1 = np.full(l-1, 1 / pow(h2, 2))
  78.         C1 = np.full(l, -(2 / tau) - (2 / pow(h2, 2)))
  79.         C1[0] = 1.
  80.         B1[0] = 0.
  81.         A1[-1] = 0
  82.         C1[-1] = 1.
  83.         u_half[i] = progonka(A1, B1, C1, right_vector2.copy()[i])
  84.     a = np.zeros((l,l))
  85.     b = np.zeros((l,l))
  86.     for i in range(1,l-1):
  87.         for j in range(1,l-1):
  88.             a[i][j] = (u_half.copy()[i-1][j]-2*u_half.copy()[i][j]+u_half.copy()[i+1][j])/pow(h1,2)
  89.             b[i][j] = (u_half.copy()[i][j-1]-2*u_half.copy()[i][j]+u_half.copy()[i][j+1])/pow(h2, 2)
  90.             u_n1[i][j] = u[i][j]+tau*(a[i][j]+b[i][j])-tau*Function.copy()[i][j]
  91.     return u_n1
  92.  
  93.  
  94.  
  95. def fill_array(length,x,y):#возвращает двумерный массив с известными границами
  96.     u = np.zeros((length,length))
  97.     for i in range(length):
  98.         u[0][i] = border(3,x[1][i],y[i,0])
  99.         u[l-1][i] = border(4,x[1][i],y[i,0])
  100.         u[i][0] = border(1,x[0][i],y[i,0])
  101.         u[i][l-1] = border(2,x[0][i],y[i,0])
  102.     return u
  103.  
  104.  
  105.  
  106. def FindEstTime_and_last_layer(epsilon,U_p,U_n):
  107.     layer_num = 0
  108.     while np.amax(np.abs(U_p-U_n)) > epsilon:
  109.         U_p = U_n
  110.         U_n = solution(U_p,X,Y,len(X),len(Y))
  111.         layer_num = layer_num + 1
  112.         #print(layer_num)
  113.     return U_n
  114.  
  115. def find_layer_num(u,u_next,time):
  116.     k = time/tau
  117.     i=1
  118.     while(k-i>0):
  119.         print(k-i)
  120.         u=u_next
  121.         u_next = solution(u.copy(),X,Y)
  122.         i=i+1
  123.     return u_next
  124.  
  125. NX=40
  126. NY=40
  127. X = np.linspace(0.,1.,NX)
  128. Y = np.linspace(0.,1.,NY)
  129. X, Y = np.meshgrid(X, Y)
  130. l = len(X)
  131. h = 1/l
  132. tau = pow(h,2)/sin(pi*h)
  133. epsilon = pow(10,-6)
  134. u = fill_array(l,X,Y)
  135. u_next = solution(u,X,Y,len(X),len(Y))
  136. u_last = FindEstTime_and_last_layer(epsilon,u.copy(),u_next.copy())
  137. Z = exactsolution(X,Y)
  138. print(np.amax(abs(u_last-Z)))
  139. fig = plt.figure(figsize=plt.figaspect(1))
  140. ax = fig.add_subplot(1, 1, 1, projection='3d')
  141. surf1 = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, cmap='Reds')
  142. surf2 = ax.plot_surface(X, Y, u_last, rstride=1, cstride=1, linewidth=0, cmap='Greens')
  143. plt.show()
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