Advertisement
HellFinger

Untitled

Mar 29th, 2022
802
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # 4-ая лаба 2-ой вариант
  2. def u_0yt(x,y,t):
  3.     return 0
  4.  
  5. def u_1yt(x,y,t):
  6.     return sin(y)*sin(t)
  7.  
  8. def u_x0t(x,y,t):
  9.     return 0
  10.  
  11. def uy_x1t(x,y,t):
  12.     return -sin(x)*sin(t)
  13.  
  14. def u_xy0(x,y,t):
  15.     return 0
  16.  
  17. def f_xyt(x,y,t):
  18.     return sin(x)*sin(y)*(cos(t)+2*sin(t))
  19.  
  20. def main3(u_0yt,u_1yt,u_x0t,u_x1t,u_xy0,x0,x1,y0,y1,h,t0,t1,tau):
  21.     In = int((x1-x0)/h)
  22.     Jn = int((y1-y0)/h)
  23.     Kn = int((t1-t0)/tau)
  24.     tau_2 = tau/2
  25.     r=tau/(2*h*h)
  26.    
  27.     I = np.linspace(x0,x1,In) # массив с x
  28.     J = np.linspace(y0,y1,Jn) # массив с y
  29.     K = np.linspace(t0,t1,Kn) # массив с t
  30.    
  31.     w=np.zeros((Jn,In,Kn))
  32.    
  33.     # 0yt уже заполнен нулями
  34.     # 1yt
  35.     for k in range(0,Kn):
  36.         for j in range(0,Jn):
  37.             w[j,In-1,k] = u_1yt(x1,J[j],K[k])
  38.            
  39.     # x0t уже заполнен нулями
  40.     # xy0 уже заполнен нулями        
  41.    
  42.     for k in range (1,Kn): # идем по t
  43.  
  44.         P=np.zeros((Jn,In)) # промежуточный слой
  45.        
  46.         for j in range (1,Jn-1): # идем по y
  47.            
  48.             A=np.zeros((In-2,In-2))
  49.             b=np.zeros(In-2)
  50.            
  51.             for n in range (0,In-2):
  52.                 A[n,n] = -(2*r+1)
  53.             for n in range (0,In-3):          
  54.                 A[n+1,n] = r
  55.                 A[n,n+1] = r
  56.             for i in range (0,In-2):
  57.                 b[i] = -r*w[j+1,i,k-1]+(2*r-1)*w[j,i,k-1]-r*w[j-1,i,k-1]-r*h*h*f_xyt(I[i],J[j],K[k]-tau_2)
  58.            
  59.             #b[0] -= r*u_0yt(x0,J[j],K[k]-tau_2) вернёт 0, поэтому бесполезное действие
  60.             b[In-3] -= r*u_1yt(x1,J[j],K[k]-tau_2)
  61.            
  62.             p = solve(A, b)
  63.             #print(p[0])
  64.             for i in range(0,In-2):
  65.                 P[j,i+1] = p[i]
  66.                
  67.             #P[j,0] = u_0yt(I[0],J[j],K[k]-tau_2)
  68.             P[j,In-1] = u_1yt(I[In-1],J[j],K[k]-tau_2)
  69.        
  70.         #краевое первое не требуется тк оно 0
  71.         '''
  72.        for i in range(In):
  73.            P[0,i] = u_x0t(i,Jn-1,K[k]-tau_2)
  74.        
  75.        '''
  76.         for i in range(In):
  77.             P[Jn-1,i] = P[Jn-2,i]+h*uy_x1t(i,Jn-1,K[k]-tau_2)
  78.         for i in range (1,In-1):   # идем по x
  79.            
  80.             A=np.zeros((Jn-1,Jn-1))
  81.             b=np.zeros(Jn-1)
  82.            
  83.             for n in range (0,Jn-1):
  84.                 A[n,n]=(1+2*r)
  85.             for n in range (0,Jn-2):          
  86.                 A[n+1,n]=-r
  87.                 A[n,n+1]=-r
  88.            
  89.             for n in range (0,Jn-1):
  90.                 A[n,Jn-3] += 1
  91.                 A[n,Jn-2] -=  -1
  92.            
  93.             for j in range (0,Jn-1):
  94.                 b[j] = r*P[j,i+1]+(1-2*r)*P[j,i]+r*P[j,i-1]+r*h*h*f_xyt(I[i],J[j],K[k]-tau_2)
  95.                
  96.             b[Jn-2] -= h*uy_x1t(I[i],J[j],K[k])
  97.            
  98.             p = solve(A, b)
  99.            
  100.             for j in range(0,Jn-1):
  101.                 w[j,i,k] = p[j]
  102.            
  103.     return I, J, K, w        
  104.            
  105. def plot_3d_time(x,y,t,W,j):
  106.     X, Y = np.meshgrid(x, y)
  107.  
  108.     fig = plt.figure(figsize=(8,6))
  109.     ax = fig.add_subplot(111, projection='3d');
  110.     # Plot a basic wireframe.
  111.     #ax.plot_wireframe(X, Y, W[:,:,j],color='b');
  112.     ls = LightSource(315, 90) #cm.gist_earth
  113.     rgb = ls.shade(W[:,:,j], cmap=cm.coolwarm, vert_exag=0.01, blend_mode='soft')
  114.     surf = ax.plot_surface(X, Y, W[:,:,j],facecolors=rgb);
  115.     #fig.colorbar(surf, shrink=0.5, aspect=10)
  116.     ax.set_title("thermal conductivity graph at time t="+str(t[j]))
  117.     ax.set_xlabel("X")
  118.     ax.set_ylabel("Y")
  119.     ax.set_zlabel("Temp")
  120.     plt.show();
  121.  
  122. X_4, Y_4, T_4, W_4 = main3(u_0yt,u_1yt,u_x0t,uy_x1t,u_xy0,0,math.pi/2,0,math.pi,0.1,0,10,0.1)
  123.  
  124. ip_w.interact(plot_3d_time,  x = ip_w.fixed(X_4),  y = ip_w.fixed(Y_4), t=ip_w.fixed(T_4), W=ip_w.fixed(W_4), j=ip_w.IntSlider(min=0, max=len(T_4)-1, step=1))
  125.  
  126. # 5-ая лаба 2-ой вариант
  127.  
  128. def u_0_y_t(x,y,t):
  129.     return cosh(y)*exp(-3*t)
  130.  
  131. def u_1_y_t(x,y,t):
  132.     return 0
  133.  
  134. def u_x_0_t(x,y,t):
  135.     return cos(2*x)*exp(-3*t)
  136.  
  137. def u_x_1_t(x,y,t):
  138.     return (5/4)*cos(2*x)*exp(-3*t)
  139.  
  140. def u_x_y_0(x,y,t):
  141.     return cos(2*x)*cosh(y)
  142.  
  143. def main4(u_0_y_t, u_1_y_t, u_x_0_t, u_x_1_t, u_x_y_0, x0, x1, y0, y1, h, t0, t1, tau):
  144.     In = int((x1-x0)/h)
  145.     Jn = int((y1-y0)/h)
  146.     Kn = int((t1-t0)/tau)
  147.     tau_2 = tau/2
  148.     r=tau/(2*h*h)
  149.    
  150.     I = np.linspace(x0,x1,In) # массив с x
  151.     J = np.linspace(y0,y1,Jn) # массив с y
  152.     K = np.linspace(t0,t1,Kn) # массив с t
  153.    
  154.     w=np.zeros((Jn,In,Kn))
  155.    
  156.     # 0yt
  157.     for k in range(Kn):
  158.         for j in range(Jn):
  159.             w[j, 0, k] = u_0_y_t(0, J[j], K[k])
  160.        
  161.     # 1yt  уже заполнен нулями            
  162.    
  163.     # x0t
  164.     for k in range(Kn):
  165.         for i in range(In):
  166.             w[0, i, k] = u_x_0_t(I[i], 0, K[k])
  167.    
  168.     # xy0
  169.     for j in range(Jn):
  170.         for i in range(In):
  171.             w[j, i, 0] = u_x_y_0(I[i], J[j], 0)
  172.            
  173.     #x1t
  174.     for k in range(Kn):
  175.         for i in range(In):
  176.             w[Jn-1, i, k] = u_x_1_t(I[i], Jn-1, K[k])
  177.    
  178.    
  179.    
  180.     for k in range (1,Kn): # идем по t
  181.  
  182.         P=np.zeros((Jn,In)) # промежуточный слой
  183.        
  184.         for j in range (1,Jn-1): # идем по y
  185.            
  186.             A=np.zeros((In-2,In-2))
  187.             b=np.zeros(In-2)
  188.            
  189.             for n in range (0,In-2):
  190.                 A[n,n]= 2*r+1
  191.             for n in range (0,In-3):          
  192.                 A[n+1,n]= -r
  193.                 A[n,n+1]= -r
  194.             for i in range (0,In-2):
  195.                 b[i] = w[j,i+1,k-1]
  196.                
  197.             b[0] += r*u_0_y_t(I[0],J[j],K[k]-tau_2)
  198.             b[In-3] += r*u_1_y_t(I[In-1],J[j],K[k]-tau_2)
  199.                      
  200.             p = solve(A, b)
  201.            
  202.             for i in range(0,In-2):
  203.                 P[j,i+1] = p[i]
  204.                
  205.             #P[j,0] = u_0_y_t(I[0],J[j],K[k]-tau_2)
  206.             #P[j,In-1] = u_1_y_t(I[In-1],J[j],K[k]-tau_2)
  207.            
  208.         for j in range (1,Jn-1):   # идем по y 2 раз
  209.            
  210.             A=np.zeros((In-2,In-2))
  211.             b=np.zeros(In-2)
  212.            
  213.             for n in range (0,In-2):
  214.                 A[n,n]=(2*r+1)
  215.             for n in range (0,In-3):          
  216.                 A[n+1,n]=-r
  217.                 A[n,n+1]=-r
  218.            
  219.             for i in range (0,In-2):
  220.                 b[i] = P[j,i+1]
  221.                
  222.             b[0] -= -r*u_0_y_t(I[0],J[j],K[k])
  223.            
  224.             b[In-3] -= -r*u_1_y_t(I[In-1],J[j],K[k])
  225.            
  226.             p = solve(A, b)
  227.            
  228.             for i in range(0,In-3):
  229.                 w[j,i+1,k] = p[i]
  230.                
  231.         """  
  232.        for i in range (1,In-1):   # идем по x
  233.            
  234.            A=np.zeros((Jn-2,Jn-2))
  235.            b=np.zeros(Jn-2)
  236.            
  237.            for n in range (1,Jn-2):
  238.                A[n,n] = 1+2*r
  239.            for n in range (0,Jn-3):          
  240.                A[n+1,n] = -r
  241.                A[n,n+1] = -r
  242.            
  243.            for j in range (0,Jn-2):
  244.                b[j] = P[j+1,i]
  245.          
  246.            b[0] += r*w[j,0,k]
  247.            b[Jn-3] += r*w[j,In-1,k]
  248.            
  249.        
  250.            p = solve(A, b)
  251.            
  252.            for j in range(0,Jn-2):
  253.                w[j+1,i,k] = p[j]
  254.           """
  255.        
  256.     return I, J, K, w        
  257.            
  258. X_5, Y_5, T_5, W_5 = main4(u_0_y_t, u_1_y_t, u_x_0_t, u_x_1_t, u_x_y_0, 0, pi/4, 0, log(2), 0.01, 0, 3, 0.1)            
  259.  
  260. ip_w.interact(plot_3d_time,  x = ip_w.fixed(X_5),  y = ip_w.fixed(Y_5), t=ip_w.fixed(T_5), W=ip_w.fixed(W_5), j=ip_w.IntSlider(min=0, max=len(T_5)-1, step=1))
  261.  
  262. # 6-ая лаба 2-ой вариант
  263.  
  264. def K(x,s):
  265.     return (x**2)*exp(x*s)
  266.  
  267. def func(x):
  268.     return 1-x*(exp(x)-exp(-x))
  269.  
  270. def Fredholm_2_Rect(K, func, a, b, h):
  271.     n = int((b-a)/h)
  272.     x = np.linspace(a,b,n)
  273.     wt = 0.5
  274.     wj = 1
  275.     A = np.zeros((n,n))
  276.     b = np.zeros(n)
  277.     for i in range(n):
  278.         A[i,0] = -h*wt*K(x[i],x[0])
  279.         for j in range(1,n-1):
  280.             A[i,j] = -h*wj*K(x[i],x[j])
  281.         A[i,n-1] = -h*wt*K(x[i],x[n-1])
  282.         A[i,i] += 1
  283.         b[i] = func(x[i])
  284.     y = solve(A, b)
  285.     return x , y
  286.  
  287. def main5(K, func, a, b, h):
  288.     x, y = Fredholm_2_Rect(K, func, a, b, h)
  289.     fig = plt.figure(figsize=(16,8))
  290.     plt.plot(x,y,label='')
  291.     plt.title('Approximate solution',fontsize=24)
  292.     plt.xlabel('x')
  293.     plt.ylabel('u(x)')
  294.     plt.show()
  295.  
  296. main5(K, func, -1, 1, 0.001)
Advertisement
RAW Paste Data Copied
Advertisement