C++: Method of Minimum Squares used for System Solving

[Const@ntine]

Homework Statement

Okay, this one is a bit big, and I'm attempting a translation, so I'm going to post a TL;DR version at the bottom just to be sure. Anyway, here goes:

While doing an experiment, we write down the values of a physical/natural size y, for various values of a physical/natural size x. Let's say n is the amount of (x,y) pairs which are the result of how many times we did this experiment. We know that these two sizes are connected through the following formula: y = f(x) = a*x + b
Sadly, it's easy to have lots of faults with such a way of writting our measurments. So, we need to find a linear relation between x & y (so essentially the values of a & b), which is as close as it can to (x_i,y_i).
So we need to find the formula f(x) = a*x + b that minimizes the average square fault E = √(1/n ** ∑_i=1ⁿ(f(x_i) - y_i)²)
This is equilavent to this system:

n*b + (∑_i=1ⁿx_i)*a = ∑_i=1ⁿy_i

(∑_i=1ⁿx_i)*b + (∑_i=1ⁿx_i²)*a = ∑_i=1ⁿx_i*y_i

You have to write a program in C++ & Fortran that get close to the linear relation of x & y, by using the pairs of (x_i,y_i). The amount of measurements (n) and the pairs themselves will be entered through the keyboard, by the user.

Use the Program from Exercise 2 to solve the system.

TLD;DR: Create a program that takes n amount of (x,y) pairs, entered through the keyboard, and solves the above equation system, with a & b being the unknown quantities.

PS: That floating n over at the ∑ is just the "limit", it's not Σ to the nth power.

Homework Equations

The Attempt at a Solution

Okay, as far as this program goes, I have this:

#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;

int main()
{
 int i, j;
 int n;
 double coor[50][2];

 cout << "How many measurments did you take ? ";
 cin >> n;

 for(i=1; i<=n; i++)
  for(j=1; j<=n; j++)
   {
    cout << "Enter your coordinates for coor(" << i << ") (" << j << "): ";
    cin >> coor[i][j];
   }
   
   return 0;
}

Yeah, it's nothing thus far, but I can't "solve" the input problem yet, so I'm stumped.

The program that solves the system, from a previous excercise:

#include <iostream>
using namespace std;

int main()
{
 double a11, a12, a21, a22;
 double b1, b2;
 double x, y, D, Dx, Dy;

 cout << "Give the values of a11, a12, a21 & a22: " << endl;
 cin >> a11 >> a12 >> a21 >> a22;
 cout << endl;

 cout << "Give the values of b1 & b2: " << endl;
 cin >> b1 >> b2;
 cout << endl;

 D = (a11*a22 - a12*a21);
 Dx = (b1*a22 - b2*a12);
 Dy = (b2*a11 - b1*a21);

if (D != 0)
{
 x = Dx/D;
 y = Dy/D;
 cout << "For those values, x & y are: " << x << " & " << y << endl;
}

if (D = 0)
{
 if (Dx = 0)
 {
  if (Dy = 0)
  {
   cout << "Then the equations/system have infinite solutions." << endl;
  }
 }
}

if (D = 0)
{
 if (Dx != 0)
 {
  if (Dy != 0)
   {
    cout << "The equations/system are unsolvable." << endl;
   }
  }
}

return 0;
}

It's a tad different, but the basic idea is the same.

Anyway, I'd greatly appreciate any kind of help!

 

[BvU]

That floating n over at the ∑ is just the "limit", it's not Σ to the nth power

Ah, you mean $$\sum_{i=1}^n$$ In PF that's written as $$\sum_{i=1}^n$$ !

If you enter e.g. n = 10 you want to go through the

  {
    cout << "Enter your coordinates for coor(" << i << ") (" << j << "): ";
    cin >> coor[i][j];
   }

loop how many times ? As presented, i runs from 1 to n and for each i, j runs from 1 to n too, so that is n² times !

 

[Const@ntine]

Ah, you mean $$\sum_{i=1}^n$$ In PF that's written as $$\sum_{i=1}^n$$ !

Thanks! I didn't know about that, so I just improvised. I'll keep it in mind for next time.

If you enter e.g. n = 10 you want to go through the

  {
    cout << "Enter your coordinates for coor(" << i << ") (" << j << "): ";
    cin >> coor[i][j];
   }

loop how many times ? As presented, i runs from 1 to n and for each i, j runs from 1 to n too, so that is n² times !

Yeah, I noticed the problem when I run it, but I'm out of ideas. What should happen, would be this:

E.G. n = 3:
Enter your coordinates for coor[1][1] = ...
Enter your coordinates for coor[1][2] = ...
Enter your coordinates for coor[2][1] = ...
Enter your coordinates for coor[2][2] = ...
Enter your coordinates for coor[3][1] = ...
Enter your coordinates for coor[3][2] = ...

But, the result is obviously not that. I figure the problem is somewhere in the j <= ... part. If I put in i, it disregards one half of the array. if I put in n, well, it creates many more "columns". If I put in 2 (since I want my array to have one spot for x, and one for y, essentially), it stops before I'm able to put in all of my "coordinates"/data for each measurement.

 

[willem2]

arrays start at 0, so you should have j=0; j<=1 in the for loop.

 

[Const@ntine]

arrays start at 0, so you should have j=0; j<=1 in the for loop.

I know, but I was under the impression that if I declared it, I could start at 1 or whatever other number. I've not really worked with arrays before, but I figured since I created room for a big array, even if the input for n was something like, say, 25, there'd be plenty of room to start from 1 and go coor[1][1] -> coo[1][2] etc. Correct me if I'm wrong though, like I said, I've had minimal hands-on with arrays.

 

[BvU]

Let's say n is the amount of (x,y) pairs

So you only want to go through n times. Easiest solution: one array for x and one array for y !

 

[Const@ntine]

So you only want to go through n times. Easiest solution: one array for x and one array for y !

Yeah, that's actually a nifty way around it. I figured I'd get some work done with multidimensional arrays to get the hang of them, but I suppose I should familiarize myself with the simple ones first.

Okay, it's midnight here, so I'm going to try the program tomorrow and get back to the thread. Thanks for the help so far everyone!

 

[Const@ntine]

I solved it on C++, but I have a problem with the Fortran Version, about halfway through.

First up, the working C++ version:

#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;

int main()
{
 int i, j;
 int n;
 double coorx[50];
 double coory[50];
 double sumx, sumy, sumxs, sumxy;
 double a, b;
 double D, Da, Db;

 cout << "How many measurments did you take ? ";
 cin >> n;
 cout << endl;

 sumx = 0;
 sumy = 0;
 sumxs = 0;
 sumxy = 0;

 for(i=1; i<=n; i++)
  {
   cout << "Enter the x coordinate for Measurement No " << i << " in coorx(" << i << "): ";
   cin >> coorx[i];
   sumx += coorx[i];
   sumxs += (coorx[i]*coorx[i]);
  }
 
  cout << endl << "The sum of x coordinates is: " << sumx << endl << "And the sum of squared coordinates x is: " << sumxs << endl << endl;
 
  for(j=1; j<=n; j++)
  {
   cout << "Enter the y coordinate for Measurement No " << j << " in coory(" << j << "): ";
   cin >> coory[j];
   sumy += coory[j];
  }
 
  cout << endl << "The sum of y coordinates is: " << sumy << endl << endl;
 
  for(i=1, j=1; i<=n, j<=n; i++, j++)
  {
   sumxy += coorx[i]*coory[j];
  }
   
  cout << "The sum of the products between coordinates x and y for each measurement is: " << sumxy << endl << endl;
 
 D = (n*sumxs - sumx*sumx);
 Db = (sumy*sumxs - sumxy*sumx);
 Da = (sumxy*n - sumy*sumx);

if (D != 0)
{
 b = Db/D;
 a = Da/D;
 cout << "For those values, a & b are: " << a << " & " << b << endl << endl;
}

if (D = 0)
{
 if (Da = 0)
 {
  if (Db = 0)
  {
   cout << "Then the equations/system have infinite solutions." << endl << endl;
  }
 }
}

if (D = 0)
{
 if (Da != 0)
 {
  if (Db != 0)
   {
    cout << "The equations/system are unsolvable." << endl << endl;
   }
  }
}

return 0;
}

And now here's my Fortran one:

PROGRAM SquarMeth

REAL, DIMENSION(50) :: coorx
REAL, DIMENSION(50) :: coory

INTEGER :: i, j
INTEGER :: n

REAL :: sumx, sumy, sumxs, sumxy
REAL :: D, Da, Db
REAL :: a, b

PRINT *, 'How many measurements did you take ?'
READ *, n
PRINT *

sumx = 0
sumy = 0
sumxs = 0
sumxy = 0

DO i=1,n
PRINT *,'Enter the x coordinate for Measurement No', i ,'in coorx(',i,'):'
READ *, coorx(i)
sumx = sumx + coorx(i)
sumxs = sumxs + (coorx(i)*coorx(i))
ENDDO

PRINT *
PRINT *,'The sum of x coordinates is:',sumx
PRINT *
PRINT *,'And the sum of squared coordinates x is:',sumxs
PRINT *
PRINT *

DO j=1,n
PRINT *,'Enter the y coordinate for Measurement No',j,'in coory(',j,'):'
READ *, coory(j)
sumy = sumy + coory(j)
ENDDO

PRINT *
PRINT *,'The sum of y coordinates is:',sumy
PRINT *
PRINT *

DO i=1,n
DO j=1,n
sumxy = sumxy + coorx(i)*coory(j)
ENDDO 
ENDDO 

PRINT *
PRINT *,'The sum of the products between coordinates x and y for each measurement is:',sumxy
PRINT *
PRINT *

D = (n*sumxs - sumx*sumx)
Db = (sumy*sumxs - sumxy*sumx)
Da = (sumxy*n - sumy*sumx)

IF (D /= 0) THEN
b = Db / D
a = Da / D
PRINT *, 'For these values, b & a are:', b,'&',a
END IF

IF (D == 0) THEN
 IF (Db == 0) THEN
  IF (Da == 0) THEN
   PRINT *, 'The equations/system have infinite solutions.'
  END IF
 END IF
END IF

IF (D == 0) THEN
 IF (Db /= 0) THEN
  IF (Da /= 0) THEN
   PRINT *, 'The equations/system are unsolvable.'
  END IF
 END IF
END IF

END PROGRAM SquarMeth

The problem arises when it comes to the program having to compute sumxy, which skews the results of the rest of the "variables" and messes up the program.

Any help is appreciated!

 

[BvU]

Same as before: you don't want $\displaystyle{\sum_i\sum_j\operatorname {coorx}(i) \operatorname {coory}(j) \ \ }$ (n² terms), but you want $\ \ \displaystyle\sum_i\operatorname {coorx}(i) \operatorname {coory}(i)\ \$ (only n terms)

PS wouldn't it be easier for the user to enter x,y pairs instead of first a whole x array and then a whole y array ? i.e. use

DO i=1,n
PRINT *,'Enter the x and y coordinates for Measurement No', i ,'):'
READ *, coorx(i), coory(i)

sumx = sumx + coorx(i)
sumxs = sumxs + (coorx(i)*coorx(i))
sumy = ...
sumys ...
sumxy ...
ENDDO

[edit]This way, the C version would become easier to read too

 

[Const@ntine]

Same as before: you don't want $\displaystyle{\sum_i\sum_j\operatorname {coorx}(i) \operatorname {coory}(j) \ \ }$ (n² terms), but you want $\ \ \displaystyle\sum_i\operatorname {coorx}(i) \operatorname {coory}(i)\ \$ (only n terms)

PS wouldn't it be easier for the user to enter x,y pairs instead of first a whole x array and then a whole y array ? i.e. use

DO i=1,n
PRINT *,'Enter the x and y coordinates for Measurement No', i ,'):'
READ *, coorx(i), coory(i)

sumx = sumx + coorx(i)
sumxs = sumxs + (coorx(i)*coorx(i))
sumy = ...
sumys ...
sumxy ...
ENDDO

[edit]This way, the C version would become easier to read too

Ah yeah, I got it (I tried it and it worked). Thanks!

PS: Yeah, it's a bit "drawn out", but I like things being more "detailed", with more cout(s) and whatnot. I know it's not recommended but considering these are home-programs that nobody really checks, I figure it's okay to do things a bit "my way" until I get the hang of it.

Again, thanks a ton for the help!

 

[Mark44]

if (D = 0)
{
 if (Da = 0)
 {
  if (Db = 0)
  {
   cout << "Then the equations/system have infinite solutions." << endl << endl;
  }
 }
}

if (D = 0)
{
 if (Da != 0)
 {
  if (Db != 0)
   {
    cout << "The equations/system are unsolvable." << endl << endl;
   }
  }
}

Your C++ code has some glaring errors that no one mentioned. Neither of the two outer if blocks above will ever execute. The expression D = 0 is not comparing the value of D with zero -- it is assigning 0 to D. To compare two things for equality, you need to use the equality operator, ==, not the assignment operator, =.

In your code, if D happens to be, say 3 before entering the first if block, if (D = 0) sets D to 0. Since 0 is considered false, control flows to the second outer if block. Since D has already been set to 0, that test expression is still considered false, so control flows to whatever statements are after the second outer if block, which in your code would be the return 0; statement.

 

[Const@ntine]

Your C++ code has some glaring errors that no one mentioned. Neither of the two outer if blocks above will ever execute. The expression D = 0 is not comparing the value of D with zero -- it is assigning 0 to D. To compare two things for equality, you need to use the equality operator, ==, not the assignment operator, =.

In your code, if D happens to be, say 3 before entering the first if block, if (D = 0) sets D to 0. Since 0 is considered false, control flows to the second outer if block. Since D has already been set to 0, that test expression is still considered false, so control flows to whatever statements are after the second outer if block, which in your code would be the return 0; statement.

Oh darn, you're right. It compiled and worked with certain values, that I didn't catch that part. Thanks a ton!

 

[Mark44]

Oh darn, you're right. It compiled and worked with certain values, that I didn't catch that part. Thanks a ton!

The compiler can catch only syntax errors -- errors for which you have incorrect language usage. It can't catch errors in your logic (semantic errors) such as the ones I pointed out.

Also, your whole if statement contains some duplication that could be eliminated.

if (D != 0)
{
    b = Db/D;
    a = Da/D;
    cout << "For those values, a & b are: " << a << " & " << b << endl << endl;
}
else     // D == 0
{
    if (Da == 0 && Db == 0)
    {
       cout << "Then the equations/system have infinite solutions." << endl << endl;
    }
    if (Da != 0 && Db != 0)
    {
       cout << "The equations/system are unsolvable." << endl << endl;
    }
}

I have merely cleaned up your logic. I haven't looked closely enough at the whole program to say that your logic takes care of all possibilities. In particular, the above control structure doesn't deal with the cases where Da == 0 and Db != 0, or when Da != 0 and Db == 0. In the else clause above, either both Da and Db are zero, or they are both nonzero.

 

[Const@ntine]

The "equation system" part was a remnant from a previous program that we were to use, so I did a copy/paste and didn't really look into it. Thanks for the cleaning up and the piece of advice!