I would like to show, how to solve a system of linear equations in PHP with 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.

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 ```
| 1.0 |
| 1.0 |
| 1.0 |
```

We are searching for a vecor ```
| 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
```

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

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

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

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

```
/**
* @param NumArray $matrix
* @param NumArray $vector
* @return array
*/
protected static function gaussianEliminationPivoting(NumArray $matrix, NumArray $vector)
{
$shape = $matrix->getShape();
for ($i = 0; $i < $shape[0]; $i++) {
// find pivo element
$max = abs($matrix->get($i, $i)->getData());
$maxIndex = $i;
for ($j = $i+1; $j < $shape[0]; $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[0]; $j++) {
$fac = -$matrix->get($j, $i)->getData()/$matrix->get($i, $i)->getData();
$matrix->set($j, $matrix->get($j)->add($matrix->get($i)->dot($fac)));
$vector->set($j, $vector->get($j)->add($vector->get($i)->dot($fac)));
}
}
return [
'M' => $matrix,
'b' => $vector
];
}
```

```
/**
* @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[0]);
for ($i = $shape[0]-1; $i >= 0; $i--) {
$sum = 0;
for ($j = $i+1; $j < $shape[0]; $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;
}
```

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
```

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.