[php] Polynomen berekenen

[riotrick]

Ik probeer uit een aantal meetwaardes een grafiek te tekenen. De grafiek heeft de vorm van een 5e orde polynoom.

Dus bv:

y = ax^5 + bx^4 + cx^3 + dx^2 + ex + f

Een voorbeeldje van de meetwaardes waaruit een grafiek kan bestaan:

code:

x-as    y-as
1002    0
925 10
780 20
650 30
530 40
400 35
300 33
240 37
0   55

Met excel kan je dit vrij eenvoudig doen door van de punten een spreidingsgrafiek te maken en vervolgens een trendlijn toe te voegen op basis van een 5e orde polynoom. De trendlijn functie van excel kan daarmee ook de betreffende constantes van de funtie benaderen. In dit voorbeeld zou de functie het volgende worden:

y = -3E-13x^5 + 1E-09x^4 - 2E-06x^3 + 0,001x^2 - 0,2426x + 55,119

De bijbehorende grafiek ziet er ongeveer als volgt uit:


Wanneer de constantes van de functie zijn bepaald kan ik de grafiek tekenen.

Wat ik nou nog mis is hoe de berekening van de die waardes verloopt. Goed in excel is het simpel, maar deze berekening wil ik graag in een php scriptje verwerken, zodat ik niet steeds excel hoef te gebruiken om de constantes te bepalen. Kan iemand me wat hints in de juiste richting geven hoe deze bereking in zijn werk gaat?

[ Voor 7% gewijzigd door riotrick op 09-08-2003 01:14 ]

[pimlie]

Met die polynoom kan ik je niet helpen, maar wij hebben hier weleens een script geschreven die een B-Spline curve berekend. Ik zou je misschien kunnen uitleggen hoe en wat, maar hier hebben ze het beter uitgelegd

Onderaan staat een C script wat redelijk makkelijk naar php te vertalen is.

[Opi]

y(0) = 55 = f geldt iig.

Verder kan je op Google de termen "polynomial, coefficients, calculation" proberen.
Een resultaat is iig deze.

[Edit]
Als je matlab hebt kan je misschien ook het commando polyfit eens bekijken.

[Edit 2]
Misschien is dit meer iets voor W&L, maar die keuze laat ik over aan de TS.

[Edit3]

offtopic:
Ik blijft bezig.

Kan je misschien het resultaat posten?

[ Voor 48% gewijzigd door Opi op 09-08-2003 12:17 ]

[PrinsEdje80]

Ga maar eens op zoek in de Numerical Recipies naar idd iets van spline of een polynomen fit... Suc6

[ekoopman]

Wat ik meestal doe met data fitten is gnuplot gebruiken, het is ietwat omslachtig, maar gnuplot heeft een van de betere fit routines, is gratis, simpel in gebruik en snel. In jouw geval zal je eerst een kort scriptje genereren, met een system() call het scriptje uitvoeren en dan de output parsen (vrij triviaal in dit geval, met tail kom je een heel end ).

Het voordeel hiervan boven zelf een fit routine schrijven is dat het zo goed als altijd werkt en je helemaal niets van fit theorie hoeft te weten

[curry684]

OpifexMaximus schreef op 09 August 2003 @ 01:25:
Misschien is dit meer iets voor W&L, maar die keuze laat ik over aan de TS.

Afaics is de essentie:

Wat ik nou nog mis is hoe de berekening van de die waardes verloopt. Goed in excel is het simpel, maar deze berekening wil ik graag in een php scriptje verwerken, zodat ik niet steeds excel hoef te gebruiken om de constantes te bepalen.

Het gaat er dus om hoe je gestandaardiseerde wiskunde basaal uit kunt splitsen tot iets wat PHP kan volgen (erg afaics dus ).

[Opi]

curry684 schreef op 09 August 2003 @ 20:37:
[...]

Afaics is de essentie:

[...]

Het gaat er dus om hoe je gestandaardiseerde wiskunde basaal uit kunt splitsen tot iets wat PHP kan volgen (erg afaics dus ).

Ik vermoed dat de TS niet zozeer moeite heeft met PHP, maar meer met de wiskundige kant. Iets waar men in W&L meer ervaring mee heeft.

[Xerras]

Het bereken van dit soort problemen wordt vaak gedaan met een beste kwadratisch fit. Zoek eens op "least square fit" er is genoeg over te vinden. Als je niet zo wiskundig bent aangelegt dan leg ik het graag op request uit

[ Voor 10% gewijzigd door Xerras op 09-08-2003 23:19 ]

[Verwijderd]

Om even terug te komen op het idee van OpifexMaximus om Matlab te gebruiken: zelfs als je Matlab niet hebt loont het zich om via de website van Matlab het algoritme te achterhalen.

Voor het commando polyfit is dit http://www.mathworks.com/...techdoc/ref/polyfit.shtml

Ook kun je nog eens zoeken op 'polyval'

[riotrick]

Ik ben er inmiddels uit. hier de code waarmee ik de constantes van elke n'de orde polynoom kan uitrekeken:

PHP:

 function calcpolynomialleastsquares($degree, $xvals, $yvals) {

    $pointcount = count($xvals);
    
    
    for ( $row = 0; $row <= $degree; $row++ )
    {
        // Find the coeffsicients for the columns.
        for ( $col = 0; $col <= $degree; $col++ )
        {
            // reset total
            $total = 0;
            
            // Find Sum(Xi^(row + col)) over all i.
            for ( $i = 0; $i < $pointcount; $i++ )
                $total += pow( $xvals[$i], $row + $col );

            $coeffs[$row][$col] = $total;
        }

        // Find the constant term.
        $total = 0;
        for ( $i = 0; $i < $pointcount; $i++ ) {
            $total += $yvals[$i] * pow( $xvals[$i], $row );
        }

        $coeffs[$row][$degree + 1] = $total;
    }


    // attempt a gaussian elimination. 
    $max_row = count($coeffs);
    $max_col = count($coeffs[0]);
    
    // index variables
    
    for ( $row = 0; $row < $max_row; $row++ )
    {
        // Make sure coeffs[row][row] != 0.
        $factor = $coeffs[$row][$row];
        
        if ( abs($factor) < 0.001 )
        {
            // Switch this row with one that is not
            // zero in position. Find this row.
            for ( $i = $row + 1; $i < $max_row; $i++ )
            {
                if ( abs( $coeffs[$i][$row] > 0.001 ))
                {
                    // Switch rows i and row.
                    for ( $j = 0; $j < $max_col; $j++ )
                    {
                        $tmp = $coeffs[$row][$j];
                        $coeffs[$row][$j] = $coeffs[$i][$j];
                        $coeffs[$i][$j]   = $tmp;
                    }
                    
                    $factor = $coeffs[$row][$row];
                }
            }
            
            // See if we found a good row.
            if ( abs($factor) < 0.001 )
                // there was no good row, and hence no solution.
                die("Crap...");
        }

        // divide each entry in this row by coeffs[row][row]
        for ( $i = 0; $i < $max_col; $i++ ) {
            $coeffs[$row][$i] /= $factor;
        }

        // subtract this row from the others.
        for ( $i = 0; $i < $max_row; $i++ )
        {
            if ( $i != $row )
            {           
    // See what factor we will multiply
                // by before subtracting for this row.
                $factor = $coeffs[$i][$row];
                
                for ( $j = 0; $j < $max_col; $j++ ) {
                    $coeffs[$i][$j] -= $factor * $coeffs[$row][$j];
                }
            }
        }
    }

    // There is a solution.
    // grab the last column of coeffs and store it
    for ( $row = 0; $row <= $degree; $row++ ) {
        $avals[$row] = $coeffs[$row][$degree + 1];
    }

    // return the results
    return $avals;


 }


$xvals[0] = 3300;
$xvals[1] = 3170;
$xvals[2] = 3000;
$xvals[3] = 2850;
$xvals[4] = 2650;
$xvals[5] = 2520;
$xvals[6] = 2250;
$xvals[7] = 2000;
$xvals[8] = 1800;
$xvals[9] = 1600;
$xvals[10] = 1420;
$xvals[11] = 1300;
$xvals[12] = 1000;
$xvals[13] = 900;

$yvals[0] = 0;
$yvals[1] = 10;
$yvals[2] = 20;
$yvals[3] = 30;
$yvals[4] = 40;
$yvals[5] = 50;
$yvals[6] = 60;
$yvals[7] = 72;
$yvals[8] = 70;
$yvals[9] = 66;
$yvals[10] = 62;
$yvals[11] = 70;
$yvals[12] = 84;
$yvals[13] = 98;

$coefficienten = calcpolynomialleastsquares(5,$xvals,$yvals);

echo "<pre>";
print_r($coefficienten);
echo "</pre>";

[ Voor 6% gewijzigd door riotrick op 13-08-2003 21:32 ]

[Opi]

dat je hem nog ff post. Ik zal binnenkort eens grondig doorspitten.

[riotrick]

hier nog een testje van een grafiek gegenereerd in PDF met behulp van deze code. De schaal is nog niet helemaal handig, maar dat is nog wat bijschaven.

http://pdfdev.amersfoort....F_grafiek_vijfde_orde.php

[Opi]

Heb je misschien heel toevallig ook nog onderzoek gedaan naar de maximale orde die nodig is om binnen een bepaalde fout te blijven? Lijkt me namelijk wel interessant.

[riotrick]

Nee helaas. Ik wist al dat deze grafieken redelijk nauwkeurig te benaderen zijn met een 5e orde polynoom. Dus daar heb ik verder geen onderzoek naar gedaan.

Ik heb nog wel 1 vraagje. Hoe kan je van een 5e orde polynoom de snijpunten met de x-as bereken (dus waar y = 0)? Net als je met een parabool de abc-formule kunt gebruiken:

-b +/- sqrt((b^2 - (4 * a * c) / 2a))

[ Voor 42% gewijzigd door riotrick op 13-08-2003 23:53 ]

[Opi]

Hmmm, komt me bekend voor; in dit topic beschouwen ze een soortgelijk probleem met een derde orde polynoom, maar wordt er wel opgemerkt dat polynomen vanaf de vijfde orde problemen opleveren. Daarom wil ik nogmaals opmerken dat het misschien handiger is om dit topic naar W&L te verschuiven of in W&L een topic over dit probleem te beginnen. Ik heb nu even geen tijd, misschien kijk ik er binnenkort nog wel naar.