SHARE

<?php

/**
* @author www.tutorialspots.com
* @base on www.numericjs.com
*/

class numeric
{
public static function LU($A) {$n = count($A);$n1 = $n - 1;$P = array();
for ($k = 0;$k < $n; ++$k)
{
$Pk =$k;
$Ak =$A[$k];$max = abs($Ak[$k]);
for ($j =$k + 1; $j <$n; ++$j) {$absAjk = abs($A[$j][$k]); if ($max < $absAjk) {$max = $absAjk;$Pk = $j; } }$P[$k] =$Pk;
if ($Pk !=$k)
{
$A[$k] = $A[$Pk];
$A[$Pk] = $Ak;$Ak = $A[$k];
}

$Akk =$Ak[$k]; for ($i = $k + 1;$i < $n; ++$i)
{
$A[$i][$k] /=$Akk;
}
for ($i =$k + 1; $i <$n; ++$i) {$Ai = &$A[$i];
for ($j =$k + 1; $j <$n1; ++$j) {$Ai[$j] -=$Ai[$k] *$Ak[$j]; ++$j;
$Ai[$j] -= $Ai[$k] * $Ak[$j];
}
if ($j ===$n1)
$Ai[$j] -= $Ai[$k] * $Ak[$j];
}
}

return array('LU' => $A, 'P' =>$P);
}

public static function LUsolve($LUP,$b)
{
$LU =$LUP['LU'];
$n = count($LU);
$x =$b;
$P =$LUP['P'];
for ($i =$n - 1; $i !== -1; --$i)
$x[$i] = $b[$i];
for ($i = 0;$i < $n; ++$i)
{
$Pi =$P[$i]; if ($P[$i] !==$i)
{
$tmp =$x[$i];$x[$i] =$x[$Pi];$x[$Pi] =$tmp;
}

$LUi =$LU[$i]; for ($j = 0; $j <$i; ++$j) {$x[$i] -=$x[$j] *$LUi[$j]; } } for ($i = $n - 1;$i >= 0; --$i) {$LUi = $LU[$i];
for ($j =$i + 1; $j <$n; ++$j) {$x[$i] -=$x[$j] *$LUi[$j]; }$x[$i] /=$LUi[$i]; } return$x;
}

public static function solve($A,$b)
{
return self::LUsolve(self::LU($A),$b);
}
}
?>

Example:

$A = array(array(1, 2, 1), array(2, 3, 2), array(3, 4, 1));$b = array(4, 7, 8);

var_dump(numeric::solve($A,$b));

Result:

array(3) {
[0]=>
float(1)
[1]=>
float(1)
[2]=>
float(1)
}