daily pastebin goal
70%
SHARE
TWEET

Untitled

a guest Feb 20th, 2019 59 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. def Dirichlet_plane(Nx, Ny, x1, y1, x2, y2):
  2.     hx = 1/ Nx
  3.     hy = 1/ Ny
  4.    
  5.     xnodes = np.linspace(x1, x2, Nx + 1)
  6.     xnodes = np.repeat(xnodes, Ny + 1)
  7.     xnodes = np.reshape(xnodes, (Nx + 1, Ny + 1))
  8.    
  9.     ynodes = np.linspace(y1, y2, Ny + 1)
  10.     ynodes = np.repeat(ynodes, Nx + 1)
  11.     ynodes = np.reshape(ynodes, (Ny + 1, Nx + 1)).T
  12.    
  13.     # Заполняем таблицу узлов
  14.     nodes = np.stack((xnodes, ynodes), axis=-1)
  15.    
  16.     u = np.zeros((Nx + 1, Ny + 1))
  17.     f = np.zeros((Nx + 1, Ny + 1))
  18.    
  19.     # Заполняем матрицу для функции f
  20.     f[:, :] = fxy(nodes[..., 0], nodes[..., 1])
  21.     f[1, 1:Ny] += gxy(nodes[0, 1:Ny][:, 0], nodes[0, 1:Ny][:, 1]) / hx / hx
  22.    
  23.     # Заполняем граничные условия для u
  24.     u[0, :] = gxy(nodes[0, :, 0], nodes[0, :, 1])
  25.     u[Nx, :] = gxy(nodes[Nx, :, 0], nodes[Nx, :, 1])
  26.     u[:, 0] = gxy(nodes[:, 0, 0], nodes[:, 0, 1])
  27.     u[:, Ny] = gxy(nodes[:, Ny, 0], nodes[:, Ny, 1])
  28.     print(u)
  29.    
  30.     # Заполняем матрицу для \phi_k(n)
  31.     phi = ifft(f[1:Nx, 1:Ny].T).T
  32.    
  33.     # Заполняем матрицу для c_k(n)
  34.     c = np.zeros((Nx - 1, Ny + 1))
  35.     c[:, 0] = ifft(f[1:Nx, 0])
  36.     c[:, Ny] = ifft(f[1:Nx, Ny])
  37.    
  38.     # Решаем систему уравнений
  39.     for k in range(1, Nx):
  40.         b = np.array(phi[k - 1])
  41.         b[0] += c[k - 1, 0] / hy / hy
  42.         b[Ny - 2] += c[k - 1, Ny] / hy / hy
  43.        
  44.         A = np.zeros((Ny - 1, Ny - 1))
  45.         flat = np.ravel(A)
  46.         flat[Ny-1::Ny] = -1 / hy / hy
  47.         flat[1::Ny] = -1 / hy / hy
  48.         flat[0::Ny] = 2 / hy / hy + 4 * (np.sin(math.pi * k * hx / 2) ** 2)
  49.         c[k - 1, 1:Ny] = np.linalg.solve(A, b)
  50.    
  51.     # Находим функцию u
  52.     u[1:Nx, 1:Ny] = ifft(c[:, 1:Ny].T).T
  53.     return u
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