impressive_i

Sweep method

Feb 23rd, 2021
791
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # created by Kirill Halo
  2. # vk.com/impressive_i
  3.  
  4. import math
  5.  
  6. # Пробные данные для уравнения A*X = B
  7. a = [[ 10.8000, 0.0475,      0, 0     ],
  8.      [  0.0321, 9.9000, 0.0523, 0     ],
  9.      [       0, 0.0369, 9.0000, 0.0570],
  10.      [       0,      0, 0.0416, 8.1000]]
  11.  
  12. b = [12.1430, 13.0897, 13.6744, 13.8972]
  13.  
  14. # Решение, которое должно получиться:
  15. # x1 = 1,118587
  16. # x2 = 1,310623
  17. # x3 = 1,503186
  18. # x4 = 1,707983
  19.  
  20. # Вывод матрицы на экран
  21. def print_arr( string, namevec, a ):
  22.     if (type(a) == int) or (type(a) == float):
  23.         print(a)
  24.     else:
  25.         print( string )
  26.         for k in range(len(a)):  
  27.             print("{}[{}] = {:8.4f}".format(namevec, k, a[k]))
  28.            
  29.  
  30. # Проверка 3х-диаг. матрицы коэффициентов на корректность
  31. def isCorrectArray(a):
  32.     n = len(a)
  33.    
  34.     for row in range(0, n):
  35.         if( len(a[row]) != n ):
  36.             print('Не соответствует размерность')
  37.             return False
  38.        
  39.     for row in range(1, n - 1):
  40.         if(abs(a[row][row]) < abs(a[row][row - 1]) + abs(a[row][row + 1])):
  41.             print('Не выполнены условия достаточности')
  42.             return False
  43.  
  44.     if (abs(a[0][0]) < abs(a[0][1]))or(abs(a[n - 1][n - 1]) < abs(a[n - 1][n - 2])):
  45.         print('Не выполнены условия достаточности')
  46.         return False
  47.        
  48.    
  49.     for row in range(0, len(a)):
  50.         if( a[row][row] == 0 ):
  51.             print('Нулевые элементы на главной диагонали')
  52.             return False
  53.     return True
  54.  
  55.  
  56.  
  57. # Процедура нахождения решения 3-х диагональной матрицы
  58. def solution(a, b):
  59.     if( not isCorrectArray(a) ):
  60.         print('Ошибка в исходных данных')
  61.         return -1
  62.  
  63.     n = len(a)
  64.     x = [0 for k in range(0, n)] # обнуление вектора решений
  65.     print('Размерность матрицы: ',n,'x',n)
  66.    
  67.     # Прямой ход
  68.     v = [0 for k in range(0, n)]
  69.     u = [0 for k in range(0, n)]
  70.     # для первой 0-й строки
  71.     v[0] = a[0][1] / (-a[0][0])
  72.     u[0] = ( - b[0]) / (-a[0][0])
  73.     for i in range(1, n - 1): # заполняем за исключением 1-й и (n-1)-й строк матрицы
  74.         v[i] = a[i][i+1] / ( -a[i][i] - a[i][i-1]*v[i-1] )
  75.         u[i] = ( a[i][i-1]*u[i-1] - b[i] ) / ( -a[i][i] - a[i][i-1]*v[i-1] )
  76.     # для последней (n-1)-й строки
  77.     v[n-1] = 0
  78.     u[n-1] = (a[n-1][n-2]*u[n-2] - b[n-1]) / (-a[n-1][n-1] - a[n-1][n-2]*v[n-2])
  79.    
  80.     print_arr('Прогоночные коэффициенты v: ','v', v)
  81.     print_arr('Прогоночные коэффициенты u: ','u', u)
  82.    
  83.     # Обратный ход
  84.     x[n-1] = u[n-1]
  85.     for i in range(n-1, 0, -1):
  86.         x[i-1] = v[i-1] * x[i] + u[i-1]
  87.        
  88.     return x    
  89.                
  90.    
  91. # MAIN - блок программмы
  92. x = solution(a, b)  # Вызываем процедуру решение
  93. print_arr('Решение: ','x', x)
RAW Paste Data