• API
• FAQ
• Tools
• Archive
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
57.     Function = function(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.
72.         #print(datetime.now() - start_time)
73.     for i in range(1,l-1):
74.         for j in range(1,l-1):
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.
Top