Free Online Tutorials

<?php

/**
 * @author www.tutorialspots.com
 * @base on www.numericjs.com
 * @copyright 2015
 */
 
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)
}