fork download
copy 
  1. program QuadraticFit;
  2.  
  3. uses
  4.   SysUtils, Math;
  5.  
  6. const
  7.   N = 25; // 5x5 точек данных
  8.  
  9. type
  10.   TMatrix = array[1..6, 1..7] of Double;
  11.   TVector = array[1..6] of Double;
  12.  
  13. procedure SolveLinearSystem(var Matrix: TMatrix; var Solution: TVector);
  14. var
  15.   i, j, k, m: Integer;
  16.   Temp: Double;
  17. begin
  18.   // Метод Гаусса с выбором главного элемента
  19.   for k := 1 to 6 do
  20.   begin
  21.     // Выбор главного элемента
  22.     m := k;
  23.     for i := k+1 to 6 do
  24.       if Abs(Matrix[i,k]) > Abs(Matrix[m,k]) then
  25.         m := i;
  26.  
  27.     // Перестановка строк
  28.     if m <> k then
  29.       for j := k to 7 do
  30.       begin
  31.         Temp := Matrix[k,j];
  32.         Matrix[k,j] := Matrix[m,j];
  33.         Matrix[m,j] := Temp;
  34.       end;
  35.  
  36.     // Исключение переменной
  37.     for i := k+1 to 6 do
  38.     begin
  39.       Temp := Matrix[i,k]/Matrix[k,k];
  40.       for j := k to 7 do
  41.         Matrix[i,j] := Matrix[i,j] - Temp*Matrix[k,j];
  42.     end;
  43.   end;
  44.  
  45.   // Обратный ход
  46.   for i := 6 downto 1 do
  47.   begin
  48.     Solution[i] := Matrix[i,7];
  49.     for j := i+1 to 6 do
  50.       Solution[i] := Solution[i] - Matrix[i,j]*Solution[j];
  51.     Solution[i] := Solution[i]/Matrix[i,i];
  52.   end;
  53. end;
  54.  
  55. var
  56.   x, y, z: array[1..N] of Double;
  57.   // Основные суммы
  58.   Sx, Sy, Sz, Sx2, Sy2, Sxy, Sxz, Syz: Double;
  59.   // Дополнительные суммы
  60.   Sx3, Sy3, Sx2y, Sxy2, Sx4, Sy4, Sx2y2, Sx3y, Sxy3, Sx2z, Sxyz, Sy2z: Double;
  61.   A, B, C, D, E, F: Double;
  62.   Matrix: TMatrix;
  63.   Solution: TVector;
  64.   i: Integer;
  65. begin
  66.   // Инициализация данных (правильный способ)
  67.   // Строка 1 (y=1)
  68.   x[1] := 0; y[1] := 1; z[1] := 8;
  69.   x[2] := 1; y[2] := 1; z[2] := 15;
  70.   x[3] := 2; y[3] := 1; z[3] := 26;
  71.   x[4] := 3; y[4] := 1; z[4] := 41;
  72.   x[5] := 4; y[5] := 1; z[5] := 60;
  73.  
  74.   // Строка 2 (y=2)
  75.   x[6] := 0; y[6] := 2; z[6] := 11;
  76.   x[7] := 1; y[7] := 2; z[7] := 19;
  77.   x[8] := 2; y[8] := 2; z[8] := 31;
  78.   x[9] := 3; y[9] := 2; z[9] := 47;
  79.   x[10] := 4; y[10] := 2; z[10] := 67;
  80.  
  81.   // Строка 3 (y=3)
  82.   x[11] := 0; y[11] := 3; z[11] := 14;
  83.   x[12] := 1; y[12] := 3; z[12] := 23;
  84.   x[13] := 2; y[13] := 3; z[13] := 36;
  85.   x[14] := 3; y[14] := 3; z[14] := 53;
  86.   x[15] := 4; y[15] := 3; z[15] := 74;
  87.  
  88.   // Строка 4 (y=4)
  89.   x[16] := 0; y[16] := 4; z[16] := 17;
  90.   x[17] := 1; y[17] := 4; z[17] := 27;
  91.   x[18] := 2; y[18] := 4; z[18] := 41;
  92.   x[19] := 3; y[19] := 4; z[19] := 59;
  93.   x[20] := 4; y[20] := 4; z[20] := 81;
  94.  
  95.   // Строка 5 (y=5)
  96.   x[21] := 0; y[21] := 5; z[21] := 20;
  97.   x[22] := 1; y[22] := 5; z[22] := 31;
  98.   x[23] := 2; y[23] := 5; z[23] := 46;
  99.   x[24] := 3; y[24] := 5; z[24] := 65;
  100.   x[25] := 4; y[25] := 5; z[25] := 88;
  101.  
  102.   // Обнуление всех сумм
  103.   Sx := 0; Sy := 0; Sz := 0;
  104.   Sx2 := 0; Sy2 := 0; Sxy := 0;
  105.   Sxz := 0; Syz := 0;
  106.   Sx3 := 0; Sy3 := 0;
  107.   Sx2y := 0; Sxy2 := 0;
  108.   Sx4 := 0; Sy4 := 0;
  109.   Sx2y2 := 0;
  110.   Sx3y := 0; Sxy3 := 0;
  111.   Sx2z := 0; Sxyz := 0; Sy2z := 0;
  112.  
  113.   // Вычисление всех необходимых сумм
  114.   for i := 1 to N do
  115.   begin
  116.     Sx := Sx + x[i];
  117.     Sy := Sy + y[i];
  118.     Sz := Sz + z[i];
  119.     Sx2 := Sx2 + x[i]*x[i];
  120.     Sy2 := Sy2 + y[i]*y[i];
  121.     Sxy := Sxy + x[i]*y[i];
  122.     Sxz := Sxz + x[i]*z[i];
  123.     Syz := Syz + y[i]*z[i];
  124.     Sx3 := Sx3 + x[i]*x[i]*x[i];
  125.     Sy3 := Sy3 + y[i]*y[i]*y[i];
  126.     Sx2y := Sx2y + x[i]*x[i]*y[i];
  127.     Sxy2 := Sxy2 + x[i]*y[i]*y[i];
  128.     Sx4 := Sx4 + x[i]*x[i]*x[i]*x[i];
  129.     Sy4 := Sy4 + y[i]*y[i]*y[i]*y[i];
  130.     Sx2y2 := Sx2y2 + x[i]*x[i]*y[i]*y[i];
  131.     Sx3y := Sx3y + x[i]*x[i]*x[i]*y[i];
  132.     Sxy3 := Sxy3 + x[i]*y[i]*y[i]*y[i];
  133.     Sx2z := Sx2z + x[i]*x[i]*z[i];
  134.     Sxyz := Sxyz + x[i]*y[i]*z[i];
  135.     Sy2z := Sy2z + y[i]*y[i]*z[i];
  136.   end;
  137.  
  138.   // Формирование системы уравнений
  139.   Matrix[1,1] := N;    Matrix[1,2] := Sx;   Matrix[1,3] := Sy;   Matrix[1,4] := Sx2;  Matrix[1,5] := Sxy;  Matrix[1,6] := Sy2;  Matrix[1,7] := Sz;
  140.   Matrix[2,1] := Sx;   Matrix[2,2] := Sx2;  Matrix[2,3] := Sxy;  Matrix[2,4] := Sx3;  Matrix[2,5] := Sx2y; Matrix[2,6] := Sxy2; Matrix[2,7] := Sxz;
  141.   Matrix[3,1] := Sy;   Matrix[3,2] := Sxy;  Matrix[3,3] := Sy2;  Matrix[3,4] := Sxy2; Matrix[3,5] := Sy3;  Matrix[3,6] := Sx2y; Matrix[3,7] := Syz;
  142.   Matrix[4,1] := Sx2;  Matrix[4,2] := Sx3;  Matrix[4,3] := Sx2y; Matrix[4,4] := Sx4;  Matrix[4,5] := Sx3y; Matrix[4,6] := Sx2y2; Matrix[4,7] := Sx2z;
  143.   Matrix[5,1] := Sxy;  Matrix[5,2] := Sx2y; Matrix[5,3] := Sy3;  Matrix[5,4] := Sx3y; Matrix[5,5] := Sx2y2; Matrix[5,6] := Sxy3; Matrix[5,7] := Sxyz;
  144.   Matrix[6,1] := Sy2;  Matrix[6,2] := Sxy2; Matrix[6,3] := Sx2y; Matrix[6,4] := Sx2y2; Matrix[6,5] := Sxy3; Matrix[6,6] := Sy4;  Matrix[6,7] := Sy2z;
  145.  
  146.   // Решение системы
  147.   SolveLinearSystem(Matrix, Solution);
  148.  
  149.   // Присвоение коэффициентов
  150.   A := Solution[1];
  151.   B := Solution[2];
  152.   C := Solution[3];
  153.   D := Solution[4];
  154.   E := Solution[5];
  155.   F := Solution[6];
  156.  
  157.   // Вывод результатов
  158.   Writeln('Коэффициенты функции F(x,y) = A + Bx + Cy + Dx² + Exy + Fy²:');
  159.   Writeln('A = ', A:0:6);
  160.   Writeln('B = ', B:0:6);
  161.   Writeln('C = ', C:0:6);
  162.   Writeln('D = ', D:0:6);
  163.   Writeln('E = ', E:0:6);
  164.   Writeln('F = ', F:0:6);
  165. end.
Success #stdin #stdout 0s 5316KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Коэффициенты функции F(x,y) = A + Bx + Cy + Dx² + Exy + Fy²:
A = 6.278463
B = 5.668266
C = 0.089275
D = 2.000000
E = 0.443911
F = 0.677610
https://ideone.com/kVsDMT
language:
Pascal (fpc 3.0.4)
created:
18 hours 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!