# Solving a linear system in PHP

I would like to show, how to solve a system of linear equations in PHP with NumPHP.

## What is NumPHP?

NumPHP is my pretty small library in pre development status. I have build it after using Numpy in the university. I have missed some features like the multiplication of matrices in PHP. Sure, PHP is a language for the web and not for science, but sometimes it can be very helpful to solve problems with the calculation of matrices. Numpy can do a lot more and NumPHP hasn't the claim to be used for science, but maybe it can be helpful to solve little problems on a webpage or to evaluate some user datas. I have called the current status pre developer, cause I just wanted to check if it would be possible and I just have implemented the most important functions for this little example. NumPHP is completely tested. Enough excuses.

## System of linear equations

I have a pretty small problem as example. We have matrix A

``````| 11.0 , 44.0 ,  1.0 |
|  0.1 ,  0.4 ,  3.0 |
|  0.0 ,  1.0 , -1.0 |
``````
and we have vector b
``````| 1.0 |
| 1.0 |
| 1.0 |
``````
We are searching for a vecor x, so that
``````| 11.0 , 44.0 ,  1.0 |   | x1 |   | 1.0 |
|  0.1 ,  0.4 ,  3.0 | * | x2 | = | 1.0 |
|  0.0 ,  1.0 , -1.0 |   | x3 |   | 1.0 |
``````
Or as linear system in the school
``````(1) 11.0*x1 + 44.0*x2 + 1.0*x3 = 1.0
(2)  0.1*x1 +  0.4*x2 + 3.0*x3 = 1.0
(3)                x2      -x3 = 1.0
``````

## Gaussian Elimination

Maybe not the best, but a comprehensible for this problem is the Gaussian elimination. We will use it to solve your problem.

### Import NumArray

``````use NumPHP\Core\NumArray;
``````

### Create a matrix

``````\$matrixA = new NumArray(
[
[11.0, 44.0, 1.0],
[0.1, 0.4, 3.0],
[0.0, 1.0, -1.0],
]
);
``````

### Create a vector

``````\$vectorB = new NumArray(
[1.0, 1.0, 1.0]
);
``````

### Gaussian elimination with pivoting

``````/**
* @param NumArray \$matrix
* @param NumArray \$vector
* @return array
*/
protected static function gaussianEliminationPivoting(NumArray \$matrix, NumArray \$vector)
{
\$shape = \$matrix->getShape();

for (\$i = 0; \$i < \$shape; \$i++) {
// find pivo element
\$max = abs(\$matrix->get(\$i, \$i)->getData());
\$maxIndex = \$i;
for (\$j = \$i+1; \$j < \$shape; \$j++) {
\$abs = abs(\$matrix->get(\$j, \$i)->getData());
if (\$abs > \$max) {
\$max = \$abs;
\$maxIndex = \$j;
}
}
// pivoting
if (\$maxIndex !== \$i) {
// change maxIndex row with i row in \$matrix
\$temp = \$matrix->get(\$i);
\$matrix->set(\$i, \$matrix->get(\$maxIndex));
\$matrix->set(\$maxIndex, \$temp);
// change maxIndex row with i row in \$vector
\$temp = \$vector->get(\$i);
\$vector->set(\$i, \$vector->get(\$maxIndex));
\$vector->set(\$maxIndex, \$temp);
}
for (\$j = \$i+1; \$j < \$shape; \$j++) {
\$fac = -\$matrix->get(\$j, \$i)->getData()/\$matrix->get(\$i, \$i)->getData();
}
}

return [
'M' => \$matrix,
'b' => \$vector
];
}
``````

### Back substitution

``````/**
* @param NumArray \$matrix
* @param NumArray \$vector
* @return NumArray
*/
protected static function backSubstitution(NumArray \$matrix, NumArray \$vector)
{
\$shape = \$matrix->getShape();
\$xVector = \NumPHP\Core\NumPHP::zeros(\$shape);

for (\$i = \$shape-1; \$i >= 0; \$i--) {
\$sum = 0;
for (\$j = \$i+1; \$j < \$shape; \$j++) {
\$sum += \$matrix->get(\$i, \$j)->dot(\$xVector->get(\$j))->getData();
}
\$xVector->set(\$i, (\$vector->get(\$i)->getData()-\$sum)/\$matrix->get(\$i, \$i)->getData());
}

return \$xVector;
}
``````

### Solution

We only have to combine these two functions and put matrix A and vector b inside.

``````\$gaussianResult = self::gaussianEliminationPivoting(\$matrixA, \$vectorB);
\$xVector = self::backSubstitution(\$gaussianResult['M'], \$gaussianResult['b']);
``````
The complete example can be found under NumPHP-examples as Symfony Command. The command has the following output
``````Matrix A:
NumArray([
[11, 44, 1],
[0.1, 0.4, 3],
[0, 1, -1]
])
Vector b:
NumArray([1, 1, 1])
Solution of A*x=b?
Vector x:
NumArray([-5.2644376899696, 1.3313069908815, 0.33130699088146])
Time for calculation: 0.0018529891967773 sec
``````

### Checking the solution with Numpy

Just to be sure, that this is the right solution, I have created a small Numpy script with the same example.

``````import time
import numpy

A = numpy.array([[11., 44., 1.], [0.1, 0.4, 3.], [0., 1., -1.]])
b = numpy.array([1., 1., 1.])

time1 = time.time()
x = numpy.linalg.solve(A, b)
timeDiff = time.time()-time1

print(x)
print "%.16f" % timeDiff
``````
And it has the following output
``````[-5.26443769  1.33130699  0.33130699]
0.0000879764556885
``````
As you can see, NumPHP can be exact, but it is incredible slow. I have detected some problems and maybe we can save some milliseconds in the future.