fork download 
  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.   // Координаты x и y остаются прежними, z пересчитывается по формуле z = 5 + 4x + 3y + 2x² + xy
  68.   // Строка 1 (y=1)
  69.   x[1] := 0; y[1] := 1; z[1] := 5 + 4*0 + 3*1 + 2*0*0 + 0*1; // z = 8
  70.   x[2] := 1; y[2] := 1; z[2] := 5 + 4*1 + 3*1 + 2*1*1 + 1*1; // z = 15
  71.   x[3] := 2; y[3] := 1; z[3] := 5 + 4*2 + 3*1 + 2*2*2 + 2*1; // z = 26
  72.   x[4] := 3; y[4] := 1; z[4] := 5 + 4*3 + 3*1 + 2*3*3 + 3*1; // z = 41
  73.   x[5] := 4; y[5] := 1; z[5] := 5 + 4*4 + 3*1 + 2*4*4 + 4*1; // z = 60
  74.  
  75.   // Строка 2 (y=2)
  76.   x[6] := 0; y[6] := 2; z[6] := 5 + 4*0 + 3*2 + 2*0*0 + 0*2; // z = 11
  77.   x[7] := 1; y[7] := 2; z[7] := 5 + 4*1 + 3*2 + 2*1*1 + 1*2; // z = 19
  78.   x[8] := 2; y[8] := 2; z[8] := 5 + 4*2 + 3*2 + 2*2*2 + 2*2; // z = 31
  79.   x[9] := 3; y[9] := 2; z[9] := 5 + 4*3 + 3*2 + 2*3*3 + 3*2; // z = 47
  80.   x[10] := 4; y[10] := 2; z[10] := 5 + 4*4 + 3*2 + 2*4*4 + 4*2; // z = 67
  81.  
  82.   // Строка 3 (y=3)
  83.   x[11] := 0; y[11] := 3; z[11] := 5 + 4*0 + 3*3 + 2*0*0 + 0*3; // z = 14
  84.   x[12] := 1; y[12] := 3; z[12] := 5 + 4*1 + 3*3 + 2*1*1 + 1*3; // z = 23
  85.   x[13] := 2; y[13] := 3; z[13] := 5 + 4*2 + 3*3 + 2*2*2 + 2*3; // z = 36
  86.   x[14] := 3; y[14] := 3; z[14] := 5 + 4*3 + 3*3 + 2*3*3 + 3*3; // z = 53
  87.   x[15] := 4; y[15] := 3; z[15] := 5 + 4*4 + 3*3 + 2*4*4 + 4*3; // z = 74
  88.  
  89.   // Строка 4 (y=4)
  90.   x[16] := 0; y[16] := 4; z[16] := 5 + 4*0 + 3*4 + 2*0*0 + 0*4; // z = 17
  91.   x[17] := 1; y[17] := 4; z[17] := 5 + 4*1 + 3*4 + 2*1*1 + 1*4; // z = 27
  92.   x[18] := 2; y[18] := 4; z[18] := 5 + 4*2 + 3*4 + 2*2*2 + 2*4; // z = 41
  93.   x[19] := 3; y[19] := 4; z[19] := 5 + 4*3 + 3*4 + 2*3*3 + 3*4; // z = 59
  94.   x[20] := 4; y[20] := 4; z[20] := 5 + 4*4 + 3*4 + 2*4*4 + 4*4; // z = 81
  95.  
  96.   // Строка 5 (y=5)
  97.   x[21] := 0; y[21] := 5; z[21] := 5 + 4*0 + 3*5 + 2*0*0 + 0*5; // z = 20
  98.   x[22] := 1; y[22] := 5; z[22] := 5 + 4*1 + 3*5 + 2*1*1 + 1*5; // z = 31
  99.   x[23] := 2; y[23] := 5; z[23] := 5 + 4*2 + 3*5 + 2*2*2 + 2*5; // z = 46
  100.   x[24] := 3; y[24] := 5; z[24] := 5 + 4*3 + 3*5 + 2*3*3 + 3*5; // z = 65
  101.   x[25] := 4; y[25] := 5; z[25] := 5 + 4*4 + 3*5 + 2*4*4 + 4*5; // z = 88
  102.  
  103.   // Обнуление всех сумм
  104.   Sx := 0; Sy := 0; Sz := 0;
  105.   Sx2 := 0; Sy2 := 0; Sxy := 0;
  106.   Sxz := 0; Syz := 0;
  107.   Sx3 := 0; Sy3 := 0;
  108.   Sx2y := 0; Sxy2 := 0;
  109.   Sx4 := 0; Sy4 := 0;
  110.   Sx2y2 := 0;
  111.   Sx3y := 0; Sxy3 := 0;
  112.   Sx2z := 0; Sxyz := 0; Sy2z := 0;
  113.  
  114.   // Вычисление всех необходимых сумм
  115.   for i := 1 to N do
  116.   begin
  117.     Sx := Sx + x[i];
  118.     Sy := Sy + y[i];
  119.     Sz := Sz + z[i];
  120.     Sx2 := Sx2 + x[i]*x[i];
  121.     Sy2 := Sy2 + y[i]*y[i];
  122.     Sxy := Sxy + x[i]*y[i];
  123.     Sxz := Sxz + x[i]*z[i];
  124.     Syz := Syz + y[i]*z[i];
  125.     Sx3 := Sx3 + x[i]*x[i]*x[i];
  126.     Sy3 := Sy3 + y[i]*y[i]*y[i];
  127.     Sx2y := Sx2y + x[i]*x[i]*y[i];
  128.     Sxy2 := Sxy2 + x[i]*y[i]*y[i];
  129.     Sx4 := Sx4 + x[i]*x[i]*x[i]*x[i];
  130.     Sy4 := Sy4 + y[i]*y[i]*y[i]*y[i];
  131.     Sx2y2 := Sx2y2 + x[i]*x[i]*y[i]*y[i];
  132.     Sx3y := Sx3y + x[i]*x[i]*x[i]*y[i];
  133.     Sxy3 := Sxy3 + x[i]*y[i]*y[i]*y[i];
  134.     Sx2z := Sx2z + x[i]*x[i]*z[i];
  135.     Sxyz := Sxyz + x[i]*y[i]*z[i];
  136.     Sy2z := Sy2z + y[i]*y[i]*z[i];
  137.   end;
  138.  
  139.   // Формирование системы уравнений
  140.   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;
  141.   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;
  142.   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;
  143.   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;
  144.   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;
  145.   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;
  146.  
  147.   // Решение системы
  148.   SolveLinearSystem(Matrix, Solution);
  149.  
  150.   // Присвоение коэффициентов
  151.   A := Solution[1];
  152.   B := Solution[2];
  153.   C := Solution[3];
  154.   D := Solution[4];
  155.   E := Solution[5];
  156.   F := Solution[6];
  157.  
  158.   // Вывод результатов
  159.   Writeln('Коэффициенты функции F(x,y) = A + Bx + Cy + Dx² + Exy + Fy²:');
  160.   Writeln('A = ', A:0:6);
  161.   Writeln('B = ', B:0:6);
  162.   Writeln('C = ', C:0:6);
  163.   Writeln('D = ', D:0:6);
  164.   Writeln('E = ', E:0:6);
  165.   Writeln('F = ', F:0:6);
  166. end.
Success #stdin #stdout 0.01s 5312KB
comments (?)
stdin 
Standard input is empty
stdout 
Коэффициенты функции 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/mE3p3o
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!