fork(2) download
copy 
  1. program SimpleIterationMethod;
  2.  
  3. uses
  4.   Math;
  5.  
  6. const
  7.   EPS = 1e-12;  // Точность вычислений
  8.   MAX_ITER = 1000;  // Максимальное количество итераций
  9.  
  10. type
  11.   // Типы для матриц и векторов
  12.   TMatrix = array[1..3, 1..3] of Double;
  13.   TVector = array[1..3] of Double;
  14.  
  15. var
  16.   A: TMatrix;    // Исходная матрица коэффициентов
  17.   B: TVector;    // Исходный вектор свободных членов
  18.   Alpha: TMatrix; // Матрица коэффициентов в нормальном виде
  19.   Beta: TVector;  // Вектор свободных членов в нормальном виде
  20.   X, X_prev: TVector; // Текущее и предыдущее приближения
  21.   NormAlpha: Double;  // Норма матрицы Alpha
  22.   k_apr: Integer;     // Априорная оценка количества итераций
  23.   k_real: Integer;    // Реальное количество итераций
  24.   converged: Boolean; // Флаг сходимости
  25.   apostEstimate: Double; // Апостериорная оценка
  26.  
  27. // 1. Процедура приведения системы к нормальному виду
  28. procedure ConvertToNormalForm(var Alpha: TMatrix; var Beta: TVector);
  29. var
  30.   i, j: Integer;
  31. begin
  32.   for i := 1 to 3 do
  33.   begin
  34.     // Вычисляем Beta[i] = B[i] / A[i,i]
  35.     Beta[i] := B[i] / A[i, i];
  36.  
  37.     // Вычисляем элементы матрицы Alpha
  38.     for j := 1 to 3 do
  39.     begin
  40.       if i = j then
  41.         Alpha[i, j] := 0  // Диагональные элементы равны 0
  42.       else
  43.         Alpha[i, j] := -A[i, j] / A[i, i];
  44.     end;
  45.   end;
  46. end;
  47.  
  48. // 2. Функция проверки условия сходимости (по норме матрицы)
  49. function CheckConvergence(Alpha: TMatrix): Boolean;
  50. var
  51.   i, j: Integer;
  52.   rowSum, maxSum: Double;
  53. begin
  54.   maxSum := 0;
  55.  
  56.   // Проверяем условие сходимости по строкам (кубическая норма)
  57.   for i := 1 to 3 do
  58.   begin
  59.     rowSum := 0;
  60.     for j := 1 to 3 do
  61.       rowSum := rowSum + Abs(Alpha[i, j]);
  62.  
  63.     if rowSum > maxSum then
  64.       maxSum := rowSum;
  65.   end;
  66.  
  67.   NormAlpha := maxSum;
  68.   CheckConvergence := (NormAlpha < 1);
  69.  
  70.   if not CheckConvergence then
  71.     Writeln('Внимание: условие сходимости не выполнено! Норма = ', NormAlpha:0:6);
  72. end;
  73.  
  74. // 3. Выбор начального приближения
  75. function GetInitialApproximation(Beta: TVector): TVector;
  76. var
  77.   i: Integer;
  78.   X0: TVector;
  79. begin
  80.   // В качестве начального приближения берем вектор Beta
  81.   for i := 1 to 3 do
  82.     X0[i] := Beta[i];
  83.  
  84.   GetInitialApproximation := X0;
  85. end;
  86.  
  87. // 4. Вычисление априорной оценки количества итераций
  88. function CalculateAPrioriEstimation(NormAlpha: Double; EPS: Double): Integer;
  89. var
  90.   k: Integer;
  91. begin
  92.   // k >= ln(eps * (1 - ||Alpha||) / ||Beta||) / ln(||Alpha||)
  93.   // Для упрощения используем приближенную формулу
  94.   if NormAlpha < 1 then
  95.     k := Ceil(Ln(EPS * (1 - NormAlpha)) / Ln(NormAlpha))
  96.   else
  97.     k := MAX_ITER;
  98.  
  99.   if k < 1 then k := 1;
  100.   if k > MAX_ITER then k := MAX_ITER;
  101.  
  102.   CalculateAPrioriEstimation := k;
  103. end;
  104.  
  105. // 5. Метод простых итераций
  106. function SimpleIteration(Alpha: TMatrix; Beta: TVector; var X: TVector; 
  107.                          var iterations: Integer): Boolean;
  108. var
  109.   i, j: Integer;
  110.   error, maxError: Double;
  111.   X_new: TVector;
  112. begin
  113.   iterations := 0;
  114.   SimpleIteration := True;
  115.  
  116.   repeat
  117.     Inc(iterations);
  118.  
  119.     // Сохраняем предыдущее приближение
  120.     X_prev := X;
  121.  
  122.     // Вычисляем новое приближение: X_new = Alpha * X + Beta
  123.     for i := 1 to 3 do
  124.     begin
  125.       X_new[i] := Beta[i];
  126.       for j := 1 to 3 do
  127.         X_new[i] := X_new[i] + Alpha[i, j] * X[j];
  128.     end;
  129.  
  130.     X := X_new;
  131.  
  132.     // Вычисляем погрешность (максимальное изменение по компонентам)
  133.     maxError := 0;
  134.     for i := 1 to 3 do
  135.     begin
  136.       error := Abs(X[i] - X_prev[i]);
  137.       if error > maxError then
  138.         maxError := error;
  139.     end;
  140.  
  141.     // Проверяем условие остановки
  142.     if iterations >= MAX_ITER then
  143.     begin
  144.       Writeln('Достигнуто максимальное количество итераций: ', MAX_ITER);
  145.       SimpleIteration := False;
  146.       Exit;
  147.     end;
  148.  
  149.   until (maxError < EPS);
  150. end;
  151.  
  152. // 6. Апостериорная оценка погрешности
  153. function AposterioriEstimate(Alpha: TMatrix; X, X_prev: TVector): Double;
  154. var
  155.   i, j: Integer;
  156.   normDiff, normAlpha, rowSum: Double;
  157. begin
  158.   // Вычисляем норму разности приближений
  159.   normDiff := 0;
  160.   for i := 1 to 3 do
  161.   begin
  162.     if Abs(X[i] - X_prev[i]) > normDiff then
  163.       normDiff := Abs(X[i] - X_prev[i]);
  164.   end;
  165.  
  166.   // Вычисляем норму матрицы Alpha
  167.   normAlpha := 0;
  168.   for i := 1 to 3 do
  169.   begin
  170.     rowSum := 0;
  171.     for j := 1 to 3 do
  172.       rowSum := rowSum + Abs(Alpha[i, j]);
  173.  
  174.     if rowSum > normAlpha then
  175.       normAlpha := rowSum;
  176.   end;
  177.  
  178.   // Апостериорная оценка: ||X^k - X^*|| <= (||Alpha|| * ||X^k - X^{k-1}||) / (1 - ||Alpha||)
  179.   AposterioriEstimate := (normAlpha * normDiff) / (1 - normAlpha);
  180. end;
  181.  
  182. // 7. Вывод результатов
  183. procedure PrintResults(X: TVector; k_apr, k_real: Integer; 
  184.                        apostEstimate: Double; converged: Boolean);
  185. var
  186.   i: Integer;
  187. begin
  188.   Writeln('========================================');
  189.   Writeln('РЕЗУЛЬТАТЫ РЕШЕНИЯ СЛАУ');
  190.   Writeln('========================================');
  191.  
  192.   if converged then
  193.   begin
  194.     Writeln('Решение найдено за ', k_real, ' итераций:');
  195.     for i := 1 to 3 do
  196.       Writeln('x[', i, '] = ', X[i]:0:15);
  197.  
  198.     Writeln;
  199.     Writeln('Априорная оценка количества итераций: ', k_apr);
  200.     Writeln('Реальное количество итераций: ', k_real);
  201.     Writeln('Апостериорная оценка погрешности: ', apostEstimate:0:15);
  202.   end
  203.   else
  204.   begin
  205.     Writeln('Решение не найдено с заданной точностью.');
  206.     Writeln('Последнее приближение:');
  207.     for i := 1 to 3 do
  208.       Writeln('x[', i, '] = ', X[i]:0:15);
  209.   end;
  210.  
  211.   Writeln('========================================');
  212. end;
  213.  
  214. // 8. Вывод матриц и векторов (для отладки)
  215. procedure PrintMatrix(M: TMatrix; name: string);
  216. var
  217.   i, j: Integer;
  218. begin
  219.   Writeln('Матрица ', name, ':');
  220.   for i := 1 to 3 do
  221.   begin
  222.     for j := 1 to 3 do
  223.       Write(M[i, j]:10:6, ' ');
  224.     Writeln;
  225.   end;
  226.   Writeln;
  227. end;
  228.  
  229. procedure PrintVector(V: TVector; name: string);
  230. var
  231.   i: Integer;
  232. begin
  233.   Writeln('Вектор ', name, ':');
  234.   for i := 1 to 3 do
  235.     Write(V[i]:10:6, ' ');
  236.   Writeln;
  237.   Writeln;
  238. end;
  239.  
  240. // Основная программа
  241. var
  242.   i, j: Integer;
  243.   residual: Double;
  244. begin
  245.   // Инициализация исходной матрицы A
  246.   A[1, 1] := 9.9;  A[1, 2] := -1.5; A[1, 3] := 2.6;
  247.   A[2, 1] := 0.4;  A[2, 2] := 13.6; A[2, 3] := -4.2;
  248.   A[3, 1] := 0.7;  A[3, 2] := 0.4;  A[3, 3] := 7.1;
  249.  
  250.   // Инициализация вектора B
  251.   B[1] := 0;
  252.   B[2] := 8.2;
  253.   B[3] := -1.3;
  254.  
  255.   Writeln('МЕТОД ПРОСТЫХ ИТЕРАЦИЙ');
  256.   Writeln('========================================');
  257.   Writeln('Исходная матрица A:');
  258.   PrintMatrix(A, 'A');
  259.   Writeln('Вектор свободных членов B:');
  260.   PrintVector(B, 'B');
  261.  
  262.   // 1. Приведение к нормальному виду
  263.   ConvertToNormalForm(Alpha, Beta);
  264.  
  265.   Writeln('После приведения к нормальному виду:');
  266.   PrintMatrix(Alpha, 'Alpha');
  267.   PrintVector(Beta, 'Beta');
  268.  
  269.   // 2. Проверка условия сходимости
  270.   if not CheckConvergence(Alpha) then
  271.   begin
  272.     Writeln('Программа продолжает работу, но сходимость не гарантирована.');
  273.   end;
  274.  
  275.   Writeln('Норма матрицы Alpha: ', NormAlpha:0:6);
  276.  
  277.   // 3. Выбор начального приближения
  278.   X := GetInitialApproximation(Beta);
  279.   Writeln('Начальное приближение:');
  280.   PrintVector(X, 'X0');
  281.  
  282.   // 4. Априорная оценка количества итераций
  283.   k_apr := CalculateAPrioriEstimation(NormAlpha, EPS);
  284.   Writeln('Априорная оценка количества итераций: ', k_apr);
  285.  
  286.   // 5. Решение методом простых итераций
  287.   converged := SimpleIteration(Alpha, Beta, X, k_real);
  288.  
  289.   if converged then
  290.   begin
  291.     // 6. Апостериорная оценка
  292.     apostEstimate := AposterioriEstimate(Alpha, X, X_prev);
  293.  
  294.     // 7. Вывод результатов
  295.     PrintResults(X, k_apr, k_real, apostEstimate, True);
  296.  
  297.     // Дополнительная проверка: подстановка решения в исходную систему
  298.     Writeln('Проверка решения (A*X - B):');
  299.     for i := 1 to 3 do
  300.     begin
  301.       residual := 0;
  302.       for j := 1 to 3 do
  303.         residual := residual + A[i, j] * X[j];
  304.       residual := residual - B[i];
  305.       Writeln('Уравнение ', i, ': невязка = ', residual:0:15);
  306.     end;
  307.   end
  308.   else
  309.   begin
  310.     PrintResults(X, k_apr, k_real, 0, False);
  311.   end;
  312.  
  313.   Writeln('Нажмите Enter для выхода...');
  314.   Readln;
  315. end.
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
МЕТОД ПРОСТЫХ ИТЕРАЦИЙ
========================================
Исходная матрица A:
Матрица A:
  9.900000  -1.500000   2.600000 
  0.400000  13.600000  -4.200000 
  0.700000   0.400000   7.100000 

Вектор свободных членов B:
Вектор B:
  0.000000   8.200000  -1.300000 

После приведения к нормальному виду:
Матрица Alpha:
  0.000000   0.151515  -0.262626 
 -0.029412   0.000000   0.308824 
 -0.098592  -0.056338   0.000000 

Вектор Beta:
  0.000000   0.602941  -0.183099 

Норма матрицы Alpha: 0.414141
Начальное приближение:
Вектор X0:
  0.000000   0.602941  -0.183099 

Априорная оценка количества итераций: 32
========================================
РЕЗУЛЬТАТЫ РЕШЕНИЯ СЛАУ
========================================
Решение найдено за 16 итераций:
x[1] = 0.139653676939303
x[2] = 0.528835522671920
x[3] = -0.226660814496702

Априорная оценка количества итераций: 32
Реальное количество итераций: 16
Апостериорная оценка погрешности: 0.000000000000337
========================================
Проверка решения (A*X - B):
Уравнение 1: невязка = -0.000000000000208
Уравнение 2: невязка = -0.000000000000014
Уравнение 3: невязка = -0.000000000000307
Нажмите Enter для выхода...
https://ideone.com/vdvRUH
language:
Pascal (fpc 3.0.4)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!