fork download
copy 
  1. program QuadraticFit;
  2.  
  3. uses
  4.   SysUtils, Math;
  5.  
  6. const
  7.   N = 25; // 5x5 data points
  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.   // Gaussian elimination with partial pivoting
  19.   for k := 1 to 6 do
  20.   begin
  21.     // Find pivot
  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.     // Swap rows
  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.     // Eliminate column
  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.   // Back substitution
  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.   Sx, Sy, Sz, Sx2, Sy2, Sxy, Sxz, Syz: Double;
  58.   Sx3, Sy3, Sx2y, Sxy2, Sx4, Sy4, Sx2y2: Double;
  59.   Sx3y, Sxy3, Sx2z, Sxyz, Sy2z: Double;
  60.   A, B, C, D, E, F: Double;
  61.   Matrix: TMatrix;
  62.   Solution: TVector;
  63.   i: Integer;
  64.  
  65. begin
  66.  
  67.   x[1] := 0; y[1] := 1; z[1] := 5 + 4*0 + 3*1 + 2*0*0 + 0*1; 
  68.   x[2] := 1; y[2] := 1; z[2] := 5 + 4*1 + 3*1 + 2*1*1 + 1*1; 
  69.   x[3] := 2; y[3] := 1; z[3] := 5 + 4*2 + 3*1 + 2*2*2 + 2*1; 
  70.   x[4] := 3; y[4] := 1; z[4] := 5 + 4*3 + 3*1 + 2*3*3 + 3*1;
  71.   x[5] := 4; y[5] := 1; z[5] := 5 + 4*4 + 3*1 + 2*4*4 + 4*1;
  72.   x[6] := 0; y[6] := 2; z[6] := 5 + 4*0 + 3*2 + 2*0*0 + 0*2;
  73.   x[7] := 1; y[7] := 2; z[7] := 5 + 4*1 + 3*2 + 2*1*1 + 1*2;
  74.   x[8] := 2; y[8] := 2; z[8] := 5 + 4*2 + 3*2 + 2*2*2 + 2*2;
  75.   x[9] := 3; y[9] := 2; z[9] := 5 + 4*3 + 3*2 + 2*3*3 + 3*2; 
  76.   x[10] := 4; y[10] := 2; z[10] := 5 + 4*4 + 3*2 + 2*4*4 + 4*2; 
  77.   x[11] := 0; y[11] := 3; z[11] := 5 + 4*0 + 3*3 + 2*0*0 + 0*3;
  78.   x[12] := 1; y[12] := 3; z[12] := 5 + 4*1 + 3*3 + 2*1*1 + 1*3;
  79.   x[13] := 2; y[13] := 3; z[13] := 5 + 4*2 + 3*3 + 2*2*2 + 2*3; 
  80.   x[14] := 3; y[14] := 3; z[14] := 5 + 4*3 + 3*3 + 2*3*3 + 3*3; 
  81.   x[15] := 4; y[15] := 3; z[15] := 5 + 4*4 + 3*3 + 2*4*4 + 4*3; 
  82.   x[16] := 0; y[16] := 4; z[16] := 5 + 4*0 + 3*4 + 2*0*0 + 0*4; 
  83.   x[17] := 1; y[17] := 4; z[17] := 5 + 4*1 + 3*4 + 2*1*1 + 1*4; 
  84.   x[18] := 2; y[18] := 4; z[18] := 5 + 4*2 + 3*4 + 2*2*2 + 2*4; 
  85.   x[19] := 3; y[19] := 4; z[19] := 5 + 4*3 + 3*4 + 2*3*3 + 3*4;
  86.   x[20] := 4; y[20] := 4; z[20] := 5 + 4*4 + 3*4 + 2*4*4 + 4*4; 
  87.   x[21] := 0; y[21] := 5; z[21] := 5 + 4*0 + 3*5 + 2*0*0 + 0*5; 
  88.   x[22] := 1; y[22] := 5; z[22] := 5 + 4*1 + 3*5 + 2*1*1 + 1*5; 
  89.   x[23] := 2; y[23] := 5; z[23] := 5 + 4*2 + 3*5 + 2*2*2 + 2*5; 
  90.   x[24] := 3; y[24] := 5; z[24] := 5 + 4*3 + 3*5 + 2*3*3 + 3*5;
  91.   x[25] := 4; y[25] := 5; z[25] := 5 + 4*4 + 3*5 + 2*4*4 + 4*5; 
  92.  
  93.   // Initialize sums
  94.   Sx := 0; Sy := 0; Sz := 0; Sx2 := 0; Sy2 := 0; Sxy := 0; Sxz := 0; Syz := 0;
  95.   Sx3 := 0; Sy3 := 0; Sx2y := 0; Sxy2 := 0; Sx4 := 0; Sy4 := 0; Sx2y2 := 0;
  96.   Sx3y := 0; Sxy3 := 0; Sx2z := 0; Sxyz := 0; Sy2z := 0;
  97.  
  98.   // Compute sums
  99.   for i := 1 to N do
  100.   begin
  101.     Sx := Sx + x[i];
  102.     Sy := Sy + y[i];
  103.     Sz := Sz + z[i];
  104.     Sx2 := Sx2 + x[i]*x[i];
  105.     Sy2 := Sy2 + y[i]*y[i];
  106.     Sxy := Sxy + x[i]*y[i];
  107.     Sxz := Sxz + x[i]*z[i];
  108.     Syz := Syz + y[i]*z[i];
  109.     Sx3 := Sx3 + x[i]*x[i]*x[i];
  110.     Sy3 := Sy3 + y[i]*y[i]*y[i];
  111.     Sx2y := Sx2y + x[i]*x[i]*y[i];
  112.     Sxy2 := Sxy2 + x[i]*y[i]*y[i];
  113.     Sx4 := Sx4 + x[i]*x[i]*x[i]*x[i];
  114.     Sy4 := Sy4 + y[i]*y[i]*y[i]*y[i];
  115.     Sx2y2 := Sx2y2 + x[i]*x[i]*y[i]*y[i];
  116.     Sx3y := Sx3y + x[i]*x[i]*x[i]*y[i];
  117.     Sxy3 := Sxy3 + x[i]*y[i]*y[i]*y[i];
  118.     Sx2z := Sx2z + x[i]*x[i]*z[i];
  119.     Sxyz := Sxyz + x[i]*y[i]*z[i];
  120.     Sy2z := Sy2z + y[i]*y[i]*z[i];
  121.   end;
  122.  
  123.  
  124.   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;
  125.   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;
  126.   Matrix[3,1] := Sy;   Matrix[3,2] := Sxy;  Matrix[3,3] := Sy2;  Matrix[3,4] := Sx2y; Matrix[3,5] := Sxy2; Matrix[3,6] := Sy3;  Matrix[3,7] := Syz;
  127.   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;
  128.   Matrix[5,1] := Sxy;  Matrix[5,2] := Sx2y; Matrix[5,3] := Sxy2; Matrix[5,4] := Sx3y; Matrix[5,5] := Sx2y2; Matrix[5,6] := Sxy3; Matrix[5,7] := Sxyz;
  129.   Matrix[6,1] := Sy2;  Matrix[6,2] := Sxy2; Matrix[6,3] := Sy3;  Matrix[6,4] := Sx2y2; Matrix[6,5] := Sxy3; Matrix[6,6] := Sy4;  Matrix[6,7] := Sy2z;
  130.  
  131.  
  132.   SolveLinearSystem(Matrix, Solution);
  133.  
  134.  
  135.   A := Solution[1];
  136.   B := Solution[2];
  137.   C := Solution[3];
  138.   D := Solution[4];
  139.   E := Solution[5];
  140.   F := Solution[6];
  141.  
  142.  
  143.   Writeln('Coefficients of F(x,y) = A + Bx + Cy + Dx^2 + Exy + Fy^2:');
  144.   Writeln('A = ', A:0:6);
  145.   Writeln('B = ', B:0:6);
  146.   Writeln('C = ', C:0:6);
  147.   Writeln('D = ', D:0:6);
  148.   Writeln('E = ', E:0:6);
  149.   Writeln('F = ', F:0:6);
  150. end.
Success #stdin #stdout 0s 5312KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Coefficients of F(x,y) = A + Bx + Cy + Dx^2 + Exy + Fy^2:
A = 5.000000
B = 4.000000
C = 3.000000
D = 2.000000
E = 1.000000
F = -0.000000
https://ideone.com/zytk1M
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!