Using Visual c++ to solve systems of linear equration

[Roblikescmptrs]

Can anyone help me set-up and solve a 3 equation set of linear simultaneous equations using Visual C++. I am working on a project at work, and have taken both C++ and C, but have not written programs in a couple of years, so I am a little out of practice. The problem is this:
(A) 3X + 5Y +7Z = 10
(B) X - 2Y + Z = 0
(C) 2X + 7Y - 4Z = -5
Solving by hand I come up with;

-2*(B) -2X + 4Y -2Z = 0
add to (C) (C); 2X + 7Y - 4Z = -5
(D); 0 + 11Y - 6Z = -5

add -3*(B) (A); 3X + 5Y + 7Z = 10
to (A) -3*(B); -3X + 6Y - 3Z = 0
(E); 0 11Y + 4Z = 10

add -1*(E) (D); 11Y - 6Z = -5
to (D) -1*(E); -11Y - 4Z = -10
(F); 0 - 10Z = -15
solving for Z = 3/2

substituting 3/2 for Z in (E)
11Y + 4*(3/2) = 10
11Y + 6 = 10
11Y = 4
Y = 4/11

substituting values for y and z in (B) to obtain x;
X - 2*(4/11) + 3/2 = 0
X - 8/11 + 3/2 = 0
X = 8/11 - 3/2
X = -17/11
I checked these answers on my calculator using reduced row echelon form--rref--and they came out the same. I was going to define X, Y, and Z to be floats. I was going to somehow use A, B, and C for the three given equation, respectively, but I don't think I can do this. Can anywone help me get going on this?

 

[bcravens]

There exist routines already for solving linear equations. I would start by examining those routines. Try finding an online version of "Numerical Methods in C". This reference has code that you can essentially copy and paste into your project (once you understand what it is doing of course).

Hope this helps.

 

[Roblikescmptrs]

Thanks for the info. Rob

 

[Mechy]

You have probably got a solution by now but what you need to do is to code the Gauss Elimination method.

I had done it before (when I was a student at the university), it is quite efficient and you can even solve a 900X900 matrix if you use pointers(It takes a while though...)

Mechy.

 

[McBugzz]

First of all, using MatLab (calling the solve routine) you can get the answer: x = -17/22, y = 4/11, z = 3/2.

Now, there are indeed many routines that do all that stuff, but there's one particularly good set of progs written back in the sixties by Forsythe, Malcolm and Moler. In their book - "Computer Methods for Mathematical Computations" all this stuff is explained. In Fortran.

However, an IT student, last year I bumped into a Cpp version and fixed a few bugs.

To procedures: DECOMP and SOLVE. DECOMP does the LU-factorisation and SOLVE solves the system for any right part.

Here's what the params mean:

DECOMP:

NDIM - the largest possible matrix size (one side)
N - the dimension of your matrix (3 in this case)
A - the matrix, defined with... see the .h file
COND - this is actually a return value - big COND means that your matrix is very bad for solving
IPVT - indicates row replacements (also a result)
WORK - just a tmp vector

SOLVE:

NDIM - same as before
N - same as before
A - pass the result A from DECOMP here
B - your right part, this array will contain the solution
IPVT - pass the result IPVT from DECOMP here

matrix.h
-------------------------------------

#ifndef _MATRIX
#define _MATRIX
// changed ndim from 25 to 6 by mcbugzz on April 8, 2003
#define ndim 6
// changed from float to double by mcbugzz on April 8, 2003
#define N_TYPE double
#define FORMAT_STR "%18.11e "
#define MATRIX(name) N_TYPE name[ndim][ndim]
#define VECTOR(name,size) N_TYPE name[size]

// SOLVE( NDIM, N, A, B, IPVT )
void solve (int n,MATRIX(a),N_TYPE * b,int * ipvt);
// DECOMP( NDIM, N, A, COND, IPVT, WORK )
void decomp(int n,MATRIX(a),N_TYPE * cond, int * ipvt,N_TYPE * work);
#endif _MATRIX

marix.cpp
------------------------------------

#include <math.h>
#include &quot;matrix.h&quot;

// casting added by mcbugzz on April 7, 2003


void solve (
int n,
MATRIX(a),
N_TYPE * b ,
int * ipvt
)

{
int kb,km1,nm1,kp1,i,k,m;
N_TYPE t ;
if (n ==1) {
b[1-1]=b[1-1] / a[1-1][1-1];
return;
}
nm1=n-1;
for (k=1;k<=nm1;k++)
{
kp1=k+1;
m=ipvt[k-1];
t=b[m-1];
b[m -1]=b[k-1];
b[k-1]=t;
for (i=kp1;i<=n;i++) b[i-1]=b[i-1] + a[i-1][k-1]*t;
}
for (kb=1 ;kb<=nm1;kb++) {
km1=n-kb;
k=km1 + 1;
b[k-1]=b[k-1]/a[k-1][k-1];
t=-b[k-1];
for (i=1;i<=km1;i++) b[i-1]=b[i-1] + a[i-1][k-1]*t;
}
b[1-1]=b[1-1] / a[1-1][1-1];
}






void decomp(
int n ,
MATRIX(a),
N_TYPE * cond ,
int * ipvt ,
N_TYPE * work
)
{
int nm1,i,j,k,kp1,kb,km1,m;
N_TYPE ek,t,anorm,ynorm,znorm;
ipvt[n-1]=1;
if (n == 1) {
*cond=1.0;
if (a[1-1][1-1] != 0.0) return;
else {
*cond=static_cast<float>(1.0e+32);
return;
}
};
nm1=n-1;
anorm=0.0;
for (j=1 ;j<= n;j++)
{
t=0.0;
for (i=1 ;i<=n;i++) t=t+static_cast<float>(fabs(a[i-1][j-1]));
if (t > anorm) anorm=t;
}
for (k=1;k<=nm1;k++)
{
kp1=k + 1;
m=k;
for (i=kp1;i<=n;i++)
if (fabs(a[i-1][k-1]) > fabs(a[m-1][k-1]))
m=i;
ipvt[k-1]=m;
if (m != k) ipvt[n-1]=-ipvt[n-1];
t=a[m-1][k-1];
a[m-1][k-1]=a[k-1][k-1];
a[k-1][k-1]=t;
if (t != 0.0) {
for (i=kp1 ;i<=n;i++) a[i-1][k-1]=-a[i-1][k-1]/t;
for (j=kp1 ;j<=n;j++) {
t=a[m-1][j-1];
a[m-1][j-1]=a[k-1][j-1];
a[k-1][j-1]=t;
if (t != 0.0) for (i=kp1 ;i<=n;i++)
a[i-1][j-1]=a[i-1][j-1] + a[i-1][k-1]*t;
}
}
}
for (k=1;k<=n;k++) {
t=0.0;
if (k != 1) {
km1=k-1;
for (i=1;i<=km1;i++)
t=t + a[i-1][k-1] * work[i-1];
}
ek=1.0;
if (t < 0.0) ek=-1.0;
if (a[k-1][k-1] == 0.0) {
*cond=static_cast<float>(1.0E+20);
return;
};
work[k-1]=-(ek+t) / a[k-1][k-1];
}
for (kb=1;kb<=nm1;kb++) {
k=n-kb;
t=0.0;
kp1=k + 1;
for (i=kp1;i<=n;i++) t=t + a[i-1][k-1] * work[k-1];
work[k-1]=t;
m=ipvt[k-1];
if (m != k) {
t=work[m-1];
work[m-1]=work[k-1];
work[k-1]=t;
}
}
ynorm=0.0;
for (i=1;i<=n;i++) ynorm=ynorm + static_cast<float>(fabs(work[i-1]));
solve(n,a,work,ipvt);
znorm=0.0;
for (i=1;i<=n;i++) znorm=znorm + static_cast<float>(fabs(work[i-1]));
*cond=anorm*znorm/ynorm;
if (*cond < 1.0) *cond=1.0;
return;
}