Rendu de formules mathematiques (MAJ)


2 mars 2015

Le paquet texgd n'étant pas aisément installable sur Debian, j'ai préféré chercher autre chose et ai choisi Tex2IM.

Toutes mes formules correspondent à des éléments IMG dont les URLs pointent vers la page tex2png.php à laquelle on passe comme paramètre la formule LaTeX que ce dernier doit générer. Ce script utilise directement Tex2IM et sauvegarde l'image dans un cache :

$Formula = str_replace(Array("\\\\", "'"), Array("\\", ""), rawurldecode($_SERVER["QUERY_STRING"]));
$PicName = md5($Formula).".png";
$CacheDir = "mathcache/";
$TmpPath = $CacheDir."tmp/";
//On ne crée l'image que si elle n'existe pas déjà dans le cache
if (!file_exists($CacheDir.$PicName))       {
    //###Il faut pouvoir passer les couleurs en paramètres : actuellement elles sont en dur dans tex2im...
    $TexGDCommand = "tex2im -o '$CacheDir$PicName' '$Formula'";
    $Ret = exec($TexGDCommand, $Output);
}
//On envoie un tag HTML contenant la référence de l'image générée.
$im = imagecreatefrompng($CacheDir.$PicName);
header('Content-Type: image/png');
imagepng($im);
imagedestroy($im);

Simulation d'emprunt


20 juin 2010

Présentation & démonstration

Etant à la recherche d'un bon outil de simulation d'emprunt bancaire et n'en trouvant pas, je me suis dit que le mieux était de faire ce dont j'avais besoin. (Expliquer les paramètres)... Voici donc le résultat (j'explique les calculs après la démo).

Capital à emprunter : Mensualité maximale : Assurance : % Annuités et coûts

Notes de calculs

Quand on parle d'un emprunt à 4%, on exprime le fait que chaque année le client doit payer au prêteur 4% de la somme restant dûe à ce moment là. Ainsi, si on note les intérêts dûs pour l'année n et le capital restant dû au début de l'année n, on a : avec t le taux de l'emprunt, ici 4%, soit 0.04. La somme totale empruntée se notera K, avec .

La plupart des prêts souscrits aujourd'hui se font à "annuités fixes", ce qui se traduit par le fait que l'emprunteur remboursera chaque année (et donc chaque mois) la même somme, ce qui est bien pratique pour se projeter quand on acquiert un bien. Ce qui est également bien pratique quand on est du côté du prêteur, on sait à l'avance ce qu'il va gagner. Qui a dit "sans rien faire" ?

On déduit de ce qui précède que la somme des intérêts payés pour une année et du capital remboursé pendant cette même année est constante. Si on note le capital remboursé durant l'année n et A l'annuité fixe : .

Ainsi :

Supposons que la formule soit valide et démontrons qu'elle reste valide pour le rang n.

Les intérêts de l'emprunt pour l'année n sont : ce qui permet d'écrire que le capital remboursé pour la même année est :

De proche en proche on peut simplifier cette expression :

CQFD

Ayant exprimé et à partir des données K, A, t et p (la durée d'emprunt), chacune de ces données peut donner lieu à une question intéressante :

  • Quelle est la durée optimale pour rembourser K, connaissant t et A ?
  • Quelle est l'annuité/mensualité connaissant K, t et la durée ?
  • Quelle est le capital empruntable connaissant A, t et la durée ?
  • Etc.

Durée optimale de remboursement

On cherche à connaître p tel que .

Pour cela on a besoin de transformer la somme en formule simple, nous allons donc faire une petite digression. Or, puisque on a affaire à une suite de raison 1+t.

Si on note la somme précédente s'écrit :
On va utiliser une petite astuce :
Pour écrire :
Et au final :

La digression étant terminée, reprenons notre calcul :

Donc :

Donc, si on connait la mensualité que peut se permettre l'emprunteur, le capital dont il a besoin et le taux en vigueur pour la banque, on peut en déduire la durée d'emprunt idéale. En utilisant la formule précédente, on peut donc faire le petit outil de calcul ci-dessous :

ParamètresDurée idéale
Capital emprunté : �
Taux usuraire : %
Mensualité idéale : �

C'est l'outil idéal pour optimiser son emprunt en faisant varier légèrement les paramètres de la banque, i.e. c'est avec ça que l'on négocie !

Calcul de l'annuité

Si on fixe le capital emprunté K, ainsi que le taux et la durée on peut calculer l'annuité et donc la mensualité d'un emprunt. C'est typiquement ce qu'imprime une banque à son client en appelant ça une "simulation".

Pour calculer cela on reprend le calcul et on écrit :

ParamètresAnnuité / mensualité
Capital emprunté : �
Taux usuraire : %
Durée en années : %

Capital empruntable

Pour terminer avec ces mini-outils, si on fixe la durée, le taux et la mensualité, on peut calculer le capital empruntable par le client à l'aide du calcul suivant (en se souvenant bien que ) :

Et pour finir :

ParamètresCapital empruntable
Durée en années : %
Taux usuraire : %
Mensualité idéale : �

Représentation graphique de fonctions et exemples d'utilisation


4 août 2008

Représentation graphique de valeurs ou séries de données

A chaque fois que j'avais besoin d'afficher une courbe j'utilisais des logiciels commerciaux dédiés ou des programmes que je faisais sur le pouce mais que j'oubliais bien vite une fois leur office rempli. Fini ! J'ai à ma disposition un script PHP relativement court (ça devrait s'améliorer) qui me permet de générer des images simplement.

<?

/*
Graph de fonctions : voir en bas du source pour des exemples d'utilisations.

Pour utiliser ce module de dessins de graphes ou fonctions la fonction la plus importante est :

MathGraph_create(<données à représenter>, <chemin du fichier image produit>, <paramètres>)
*/

global  $MathGraph_fontPath;        //Le chemin vers le fichier TTF utilisé pour l'affichage des textes sur un graph.
$MathGraph_fontPath = "";

//Convertit des coordonnées logiques (x,y) du repère des fonctions à représenter en coordonnées image/pixel
function    MathGraph_funcToDisplay($x, $y, &$displayX, &$displayY, &$params)
{
    $displayX = $params["LeftMargin"]+($x-$params["LogicalXMin"])*($params["DisplayWidth"]/$params["LogicalWidth"]);
    $displayY = $params["TopMargin"] + $params["DisplayHeight"] - ($y-$params["LogicalYMin"])*($params["DisplayHeight"]/$params["LogicalHeight"]);
}

function    MathGraph_logicalMinMax(&$ValSet, &$XMin, &$XMax, &$YMin, &$YMax)
{
    $XMin = min($LogicalXs = array_keys($ValSet));
    $YMin = min($ValSet);
    $XMax = max($LogicalXs);
    $YMax = max($ValSet);
}

//Nettoie le nombre contenu dans $txt : vire les 0 à la fin de la partie décimale de la virgule
function    MathGraph_cleanNumber($txt)
{
    if (!strlen($txt))      return false;
    $DotPos = strpos($txt, ".");
    if ($DotPos === false)      return $txt;
    //
    $Limit = strlen($txt) - 1;
    while ($txt[$Limit] == "0" && $Limit > $DotPos)     $Limit--;
    if ($txt[$Limit] != ".")        $Limit++;
    $txt = substr($txt, 0, $Limit);
    return $txt;
}

//displayModes peut être une suite de mode (xcenter, ycenter etc) séparés par des '|'. Bien sûr il faut que la combinaison ait un sens.
function    MathGraph_displayText($View, $x, $y, $color, $txt, $displayModes="")
{
global  $MathGraph_fontPath;

    $DisplayX = $x;
    $DisplayY = $y;
    $BBox = imagettfbbox(12, 0, $MathGraph_fontPath, $txt);
    $Modes = explode("|", $displayModes);
    foreach($Modes as $Mode)        {
        switch ($Mode)      {
            //On calcule DisplayX de manière à ce que le centre du texte soit en $x
            case    "xcenter":
                $Width = $BBox[2] - $BBox[0];
                $DisplayX = $DisplayX - $Width/2;
                break;
            //On calcule DisplayX de manière à ce qu'il n'aille pas plus loin que $x (aligné à droite)
            case    "xright":
                $Width = $BBox[2] - $BBox[0];
                $DisplayX = $DisplayX - $Width;
                break;
            case    "ycenter":
                $Height = $BBox[1] - $BBox[7];
                $DisplayY = $DisplayY + $Height/2;
                break;
            default:    break;
        }
    }
    //
    imagettftext($View, 12, 0, $DisplayX, $DisplayY, $color, $MathGraph_fontPath, $txt);
}

//Affiche un nombre sur l'image mais en l'ayant nettoyé au préalable (enlevé les 0 à la fin qui font "sales")
function    MathGraph_displayNumber($View, $x, $y, $color, $txt, $displayMode="")
{   return MathGraph_displayText($View, $x, $y, $color, MathGraph_cleanNumber($txt), $displayMode);     }

/*Taille à représenter : la largeur pour les X ou la hauteur pour les Y, peu importe
On va déterminer un espacement qui permet segmenter la taille en au moins 2 séparations.
Cet espacement est intuitif pour les humains (ceux qui raisonnent en base 10), il ne sera jamais de 7 ou de 3, mais plutôt de 1, 2 ou 5.*/
function    MathGraph_mediumStep($size, $threshold=5)
{
    if (!$size)     $size = 1;
    //
    $SizeLog = floor(log($size, 10));
    $CurPow10 = pow(exp($SizeLog), log(10));            //Puissance de 10 courante (130=>100, 1517=>1000 etc)
    $PrevPow10 = pow(exp($SizeLog-1), log(10));         //Puissance de 10 immédiatement précédente (130=>10, 1517=>100)
    //
    $TryOuts = array($PrevPow10, $PrevPow10*2, $PrevPow10*5, $CurPow10, $CurPow10*2, $CurPow10*5);
    $PrevValue = null;
    foreach($TryOuts as $Value)     {
        //On peut mettre 3 pour les lignes grossières et 5 pour les lignes fines
        if ($size/$Value < $threshold)      return $PrevValue;
        $PrevValue = $Value;
    }
    die("MathGraph_mediumStep : je ne devrais pas être ici, mon algo est faux !");
}

function    Graph_displayValue($value)
{
    $DisplayValue = sprintf("%.08f", $value);
    if ($DisplayValue == "-0")      $DisplayValue = "0";
    return $DisplayValue;
}

//Affiche une grille indiquant les unités X et Y d'un graphique
function    MathGraph_showGrid($View, $xStep, $yStep, $color, &$params, $displayCoords=true)
{
    $XMin = $params["LogicalXMin"];
    $XMax = $params["LogicalXMax"];
    $StartX = $XMin + $xStep * (ceil(($XMax - $XMin)/2)/$xStep);
    $YMin = $params["LogicalYMin"];
    $YMax = $params["LogicalYMax"];
    $StartY = $YMin + $yStep * (ceil(($YMax - $YMin)/2)/$yStep);
    //
    $XLeft = $StartX;       //-$xStep;
    $XRight = $StartX + $xStep;
    $YBottom = $StartY;     //-$yStep;
    $YTop = $StartY + $yStep;
    while (true)        {
        $Clip = 0;
        $TextYPos = $params["height"]-$params["BottomMargin"] + 14;
        $TextXPos = $params["LeftMargin"] - 14;
        MathGraph_funcToDisplay($XLeft, $YTop, $DisplayX, $DisplayY, $params);
        if ($DisplayX < $params["width"]-$params["RightMargin"] && $DisplayX>$params["LeftMargin"])     {
            ImageLine($View, $DisplayX, $params["TopMargin"], $DisplayX, $params["height"]-$params["BottomMargin"], $color);
            if ($displayCoords)     MathGraph_displayNumber($View, $DisplayX, $TextYPos, $StrongGridColor, Graph_displayValue($XLeft), "xcenter");
        }   else    $Clip++;
        if ($DisplayY < $params["height"]-$params["BottomMargin"] && $DisplayY>$params["TopMargin"])    {
            ImageLine($View, $params["LeftMargin"], $DisplayY, $params["width"]-$params["RightMargin"], $DisplayY, $color);
            if ($displayCoords)     MathGraph_displayNumber($View, $TextXPos, $DisplayY, $StrongGridColor, Graph_displayValue($YTop), "xright|ycenter");
        }   else    $Clip++;
        MathGraph_funcToDisplay($XRight, $YBottom, $DisplayX, $DisplayY, $params);
        if ($DisplayX < $params["width"]-$params["RightMargin"] && $DisplayX>$params["LeftMargin"])     {
            ImageLine($View, $DisplayX, $params["TopMargin"], $DisplayX, $params["height"]-$params["BottomMargin"], $color);
            if ($displayCoords)     MathGraph_displayNumber($View, $DisplayX, $TextYPos, $StrongGridColor, Graph_displayValue($XRight), "xcenter");
        }   else    $Clip++;
        if ($DisplayY < $params["height"]-$params["BottomMargin"] && $DisplayY>$params["TopMargin"])    {
            ImageLine($View, $params["LeftMargin"], $DisplayY, $params["width"]-$params["RightMargin"], $DisplayY, $color);
            if ($displayCoords)     MathGraph_displayNumber($View, $TextXPos, $DisplayY, $StrongGridColor, Graph_displayValue($YBottom), "xright|ycenter");
        }   else    $Clip++;
        //
        if ($Clip == 4)     break;  //Si on n'a affiché aucune unité c'est qu'on est sorti du repère, ça ne sert à rien de continuer.
        $XLeft -= $xStep;
        $XRight += $xStep;
        $YTop += $yStep;
        $YBottom -= $yStep;
    }
}

//Transforme une couleur de type "FF00FF" en un tableau array(255, 0, 255).
function    MathGraph_parseColor($strColor)
{
    if (strlen($strColor)!=6)       return false;
    $Red = hexdec(substr($strColor, 0, 2));
    $Green = hexdec(substr($strColor, 2, 2));
    $Blue = hexdec(substr($strColor, 4, 2));
    return array($Red, $Green, $Blue);
}

function    MathGraph_time()
{
    list($usec, $sec) = explode(" ", microtime(true));
    return $sec+$usec;
}

//Calcul des largeur et hauteur logiques
function    Graph_computeExtents(&$params)
{
    $params["LogicalWidth"] = $params["LogicalXMax"] - $params["LogicalXMin"];
    if (!$params["LogicalWidth"])       $params["LogicalWidth"] = 1;
    $params["LogicalHeight"] = $params["LogicalYMax"] - $params["LogicalYMin"];
    if (!$params["LogicalHeight"])      $params["LogicalHeight"] = 1;
}

//TODO : Documenter le tableau params que l'on peut passer à la fonction pour customiser l'apparence du rendu produit.
//$func est un tableau de tableaux, où chaque sous-tableau contient les données à représenter et ce sous deux formes possibles :
//- forme simple : les sous-tableaux contiennent les données uniquement.
//- forme avancée : les sous-tableaux contiennent des paramètres à des indices prédéfinis (pour l'instant "color" et "list") dont les données à représenter à l'indice "list".
function    MathGraph_create($func, $picPath, $params=null)
{
    $Start = MathGraph_time();
    if (!strlen($picPath))      return false;
    //Inits
    if ($params == null)        $params = array();
    if (!isset($params["width"]))       $params["width"] = 500;
    $Width = $params["width"];
    if (!isset($params["height"]))      $params["height"] = 300;
    $Height = $params["height"];
    if (!isset($params["background-color"]))        $params["background-color"] = array(255,255,255);
    if (!isset($params["color"]))                   $params["color"] = array(0,0,255);
    if (!isset($params["origAxisColor"]))           $params["origAxisColor"] = array(255, 0, 0);
    if (!isset($params["origAxisWidth"]))           $params["origAxisWidth"] = 2;
global  $MathGraph_fontPath;
    $MathGraph_fontPath = isset($params["fontPath"]) ? $params["fontPath"] : "fonts/trebucbd.ttf";
    //
    if (!isset($params["LeftMargin"]))      $params["LeftMargin"] = 20;
    if (!isset($params["RightMargin"]))     $params["RightMargin"] = 20;
    if (!isset($params["TopMargin"]))       $params["TopMargin"] = 20;
    if (!isset($params["BottomMargin"]))    $params["BottomMargin"] = 20;
    //###TODO : faire des vérifications de cohérence sur les marges et la hauteur ou la largeur...
    $params["DisplayWidth"] = $Width - $params["LeftMargin"] - $params["RightMargin"];
    $params["DisplayHeight"] = $Height - $params["TopMargin"] - $params["BottomMargin"];
    //On alloue une image de la taille demandée
    $View = ImageCreateTrueColor($Width, $Height);
    //On alloue les couleurs qui vont être utilisées dans la suite du code
    $BackgroundColor = ImageColorAllocate($View, $params["background-color"][0], $params["background-color"][1], $params["background-color"][2]);
    ImageFilledRectangle($View, 0, 0, $Width-1, $Height-1, $BackgroundColor);
    $ForegroundColor = ImageColorAllocate($View, $params["color"][0], $params["color"][1], $params["color"][2]);
    $LightGridColor = ImageColorAllocate($View, 215, 215, 215);
    $MediumGridColor = ImageColorAllocate($View, 150, 150, 150);
    $StrongGridColor = ImageColorAllocate($View, 0, 0, 0);
    $Black = ImageColorAllocate($View, 0, 0, 0);
    $OrigAxesColor = ImageColorAllocate($View, $params["origAxisColor"][0], $params["origAxisColor"][1], $params["origAxisColor"][2]);
    //Détermination des minimas et maximas de la fonction à représenter
    for($n=0; $n<count($func); $n++)        {
        $CurFunc =& $func[$n];
        if (isset($CurFunc["list"]))        $ValSet =& $CurFunc["list"];
        else        $ValSet =& $CurFunc;
        MathGraph_logicalMinMax($ValSet, $XMin, $XMax, $YMin, $YMax);
        $params["LogicalXMin"] = $n ? min($params["LogicalXMin"], $XMin) : $XMin;
        $params["LogicalYMin"] = $n ? min($params["LogicalYMin"], $YMin) : $YMin;
        $params["LogicalXMax"] = $n ? max($params["LogicalXMax"], $XMax) : $XMax;
        $params["LogicalYMax"] = $n ? max($params["LogicalYMax"], $YMax) : $YMax;
    }
    //
    Graph_computeExtents($params);
    //
    if (isset($params["title"]))        MathGraph_displayText($View, $Width / 2, 20, $Black, $params["title"], "xcenter");
    //Les axes "intermédiaires"
    $MediumXStep = MathGraph_mediumStep($params["LogicalWidth"]);
    $MediumYStep = MathGraph_mediumStep($params["LogicalHeight"]);
    $LightXStep = MathGraph_mediumStep($MediumXStep, 3);
    $LightYStep = MathGraph_mediumStep($MediumYStep, 3);
    //Ajout d'une tolérance, d'une marge aux limites du repère logique
    $params["LogicalXMax"] += 2*$LightXStep;
    $params["LogicalXMin"] -= 2*$LightXStep;
    $params["LogicalYMax"] += 2*$LightYStep;
    $params["LogicalYMin"] -= 2*$LightYStep;
    Graph_computeExtents($params);      //Après avoir modifié les limites logiques du repère on doit recalculer son étendue.
    //Dessine la grille du repère
    MathGraph_showGrid($View, $LightXStep, $LightYStep, $LightGridColor, $params, false);
    MathGraph_showGrid($View, $MediumXStep, $MediumYStep, $MediumGridColor, $params, true);
    //Les deux axes passant par (0,0)
    MathGraph_funcToDisplay(0, 0, $X0, $Y0, $params);
    ImageSetThickness($View, $params["origAxisWidth"]);
    $DisplayVerticalOrigin = $params["LogicalXMin"] * $params["LogicalXMax"] < 0;       //Si la multiplication donne un résultat négatif c'est que l'axe des origines est visible.
    if ($DisplayVerticalOrigin)     ImageLine($View, $X0, $params["TopMargin"], $X0, $params["height"]-$params["BottomMargin"], $OrigAxesColor);
    $DisplayHorizontalOrigin = $params["LogicalYMin"] * $params["LogicalYMax"] < 0;     //Ibid.
    if ($DisplayHorizontalOrigin)   ImageLine($View, $params["LeftMargin"], $Y0, $params["width"]-$params["RightMargin"], $Y0, $OrigAxesColor);
    ImageSetThickness($View, 1);
    //Dessine les fonctions une par une sur la grille précédente
    for($n=0; $n<count($func); $n++)        {
        $CurFunc =& $func[$n];
        if (isset($CurFunc["list"]))        $ValSet =& $CurFunc["list"];
        else        $ValSet =& $CurFunc;
        //
        if (isset($CurFunc["color"]))       {
            $ColorComps = MathGraph_parseColor($CurFunc["color"]);
            $LocalForegroundColor = ImageColorAllocate($View, $ColorComps[0], $ColorComps[1], $ColorComps[2]);
        }   else    $LocalForegroundColor = $ForegroundColor;
        //Parcours toutes les valeurs de la fonction courante
        if (isset($CurFunc["width"]))       ImageSetThickness($View, $CurFunc["width"]);
        $PrevX = $PrevY = false;
        foreach($ValSet as $x=>$y)      {
            MathGraph_funcToDisplay($x, $y, $CurX, $CurY, $params);
            ImageLine($View, $PrevX === false ? $CurX : $PrevX, $PrevY === false ? $CurY : $PrevY, $CurX, $CurY, $LocalForegroundColor);
            if ($PrevX !== false)       ImageFilledArc($View, $PrevX, $PrevY, 6, 6, 0, 360, $MediumGridColor, 0);
            $PrevX = $CurX;         $PrevY = $CurY;
        }
        if ($PrevX !== false)       ImageFilledArc($View, $PrevX, $PrevY, 6, 6, 0, 360, $MediumGridColor, 0);
        ImageSetThickness($View, 1);
    }
    //On peut afficher le temps de génération de l'image (en excluant celui de la génération du fichier proprement dit)
    if ($params["showCPUTime"])     {
        $Duration = MathGraph_time() - $Start;
        MathGraph_displayText($View, 0, 15, $MediumGridColor, sprintf("%f", $Duration), "");
    }
    //On crée le fichier image et on détruit la ressource run-time qui la représente
    $PicExt = substr($picPath, strrpos($picPath, ".")+1);
    switch ($PicExt)        {
        case    "jpg":      $Res = imagejpeg($View, $picPath);  break;
        case    "png":      $Res = imagepng($View, $picPath);   break;
        default:            echo "Extension <b>$PicExt</b> inconnue.<br/>\n";   break;
    }
    imagedestroy($View);
    return true;
}

?>
Le source PHP du module destiné à générer des images à partir de données séquentielles.

Exemples d'utilisation

Sinus

Sinus+Cosinus

Données générées

Notes de lecture : champs turbulents pour les phénomènes gazeux


10 juillet 2008

Notes de lecture de l'artice Turbulent Wind Fields for Gaseous Phenomena par Joe Stam et Eugene Fiume.

On décrit un champ vectoriel issu de deux composantes :

  • grande échelle (déterministe, descriptif, prévisible)
  • échelle fine (stochastique + fourier en 4D) qui contient la turbulence recherchée.

Sommes de nombre


1 décembre 2007

Le problème

Le but de cette page est d'arriver à exprimer par une formule simple, i.e. sans sommation, l'expression suivante pour k et n donnés :

On connait tous la réponse pour les cas simples k=1 et k=2 :


On peut se demander "Pourquoi ça t'intéresse ?". Ma réponse est très simple : "parce que". C'est quelque chose qui m'attire personnellement et que je n'explique pas. J'ai appris, bien après avoir trouvé comment obtenir la formule pour n'importe quel k, que le travail avait déjà été fait de la même façon par A.W.F. Edwards au vingtième siècle et d'une autre complètement différente par Jacobi en 1631.

�a peut sembler très naïf mais je ne savais pas que cette formule existait depuis un bon moment et s'appelait la formule de Faulhaber, donc il faut considérer le texte ci-dessous comme un voyage d'un aveugle dans un labyrinthe imaginaire...

Pour citer mes sources, c'est Knuth qui m'a appris l'existence de Faulhaber.

La méthode

On veut donc calculer la somme S des N premiers entiers chacun élevé à la puissance k mais ça n'est pas évident à priori. L'idée est de passer par une intégrale, c'est toujours plus simple à calculer qu'une somme discrète, mais ça nous oblige à trouver une mystérieuse fonction f telle que :

Pour trouver f(x) et pour expliquer pourquoi les bornes de l'intégration ne sont pas celles de la somme, on établit l'équation suivante :

On aboutit, en posant , un polynôme en x de degré k, à une résolution qui fait intervenir le développement du binôme, par celui de , et des inversions de matrices triangulaires LU.

Résolution

L'équation {2} fait intervenir les matrices naturellement : quand on écrit f(x) sous la forme d'un polynôme, l'intégration, longue pour des valeurs de k élevée, mais très facile, donne Q(i+1)-Q(i), où Q(x) est la primitive de f(x), et cette différence contient des termes en (i+1)^j pour j <= k+1.

Cette différence, et donc l'intégrale, est une expression contenant les coefficients du polynôme f recherché. Cette expression devant être égale à , on sait tout de suite que le terme de plus haut degré de f est égal à 1. On a donc :

avec w un polynôme de degré k-1

Pour poursuivre la résolution on peut écrire l'équation {2} sous forme matricielle : A c = (1 0 0 ... 0), avec A la matrice contenant les coefficients du triangle de Pascal issus des développements binomiaux précédent, a le vecteur contenant les coefficients du polynôme f recherché et le membre de droite figurant l'égalité à .

Pour résoudre cette équation on va écrire la matrice A sous la forme A = LU, de manière à pouvoir obtenir :

où U et L sont aisément inversables puisques triangulaires. Je ne m'étends d'ailleurs pas ici sur cette inversion dont on peut trouver la description facilement.

Implémentation

On veut donc convertir la méthode ci-dessus en C++ mais il faut bien faire attention aux problèmes de précisions sur les nombres flottants qui peuvent surgir quand on fait des inversions de matrices ou des multiplications de polynômes successives. Ma méthode pour éviter ça ? Utiliser une classe représentant des fractions et manipuler des entiers ! Les calculs deviennent exacts (ce qui est le but) mais effectivement un peu plus lent (ce qui pour une fois ne m'embète pas du tout !).

Pour le code ci-dessous, les classes utilisées (fraction, polynôme, matrice et vecteur) appartiennent à ma librairie mathématiques en cours de construction.

void    Faulhaber(int k)
{
    //On remplit une matrice avec les coefficients du triangle de pascal
    QiMatrix<QiFraction<long long int>> Pascal(k+1);
    Pascal.null();
    QiPolynomial< QiFraction<long long int> >   Basic(k+1);
    Basic.setRaw(1, 1.0f, 1.0f);    //x+1
    QiPolynomial< QiFraction<long long int> >   Current(Basic);
    int i, j;
    for(i=0; i<=k; i++)     {
        for(j=0; j<=i; j++)     Pascal.set(k-j, k-i, Current.get(j));
        Current *= Basic;
    }
    Pascal.inverse();
    //On construit la matrice diagonale T avec k, k-1, ..., 2, 1 (qui est l'inverse de la matrice 1/k, 1/(k-1), ..., 1)
    QiMatrix< QiFraction<long long int> >   Diag(k+1);
    Diag.null();
    for(i=0; i<=k; i++)     Diag.set(i, i, (long long int)k+1-i);
    Diag *= Pascal;         //On fait L^-1 * T^-1
    //On multiplie par (1, 0, ..., 0) pour récupérer la première colonne
    QiVector< QiFraction<long long int> >   BaseVector(k+1), Result(k+1);
    BaseVector.base(0);
    Result = Diag * BaseVector;
    //On transforme ce vecteur en polynôme
    QiPolynomial< QiFraction<long long int> >   Formula(k+1);
    Formula.null();
    for(i=0; i<=k; i++)     Formula.set(k-i, Result.get(i));
    //On intègre ce polynôme entre 0 et N+1 en calculant sa primitive et en générant un deuxième polynôme dans lequel on remplace x par n+1
    Formula.primitive();
    QiPolynomial< QiFraction<long long int> >   NPlus1(1);
    NPlus1.setRaw(1, 1.0f, 1.0f);
    Formula.substitute(NPlus1);
    //Présentation des résultats
    printf("Formule finale :\n");   Formula.dumpObjects();
    printf("On divise par 1/%d :\n", k+1);
    Formula /= QiFraction<long long int>(1,k+1);
    Formula.dumpObjects();
}

La sortie de ce programme pour quelques valeurs de k ressemble à ce qui suit :

Formule finale : 1X+1
On divise par 1/1 : 1X+1
Résultats pour k=1 :
Formule finale : 1/2X^2+1/2X
On divise par 1/2 : 1X^2+1X
Résultats pour k=2 :
Formule finale : 1/3X^3+1/2X^2+1/6X
On divise par 1/3 : 1X^3+3/2X^2+1/2X
Résultats pour k=3 :
Formule finale : 1/4X^4+1/2X^3+1/4X^2
On divise par 1/4 : 1X^4+2X^3+1X^2
Résultats pour k=4 :
Formule finale : 1/5X^5+1/2X^4+1/3X^3-1/30X
On divise par 1/5 : 1X^5+5/2X^4+5/3X^3-1/6X
Résultats pour k=5 :
Formule finale : 1/6X^6+1/2X^5+5/12X^4-1/12X^2
On divise par 1/6 : 1X^6+3X^5+5/2X^4-1/2X^2
Résultats pour k=6 :
Formule finale : 1/7X^7+1/2X^6+1/2X^5-1/6X^3+1/42X
On divise par 1/7 : 1X^7+7/2X^6+7/2X^5-7/6X^3+1/6X
Résultats pour k=7 :
Formule finale : 1/8X^8+1/2X^7+7/12X^6-7/24X^4+1/12X^2
On divise par 1/8 : 1X^8+4X^7+14/3X^6-7/3X^4+2/3X^2
Résultats pour k=8 :
Formule finale : 1/9X^9+1/2X^8+2/3X^7-7/15X^5+2/9X^3-1/30X
On divise par 1/9 : 1X^9+9/2X^8+6X^7-21/5X^5+2X^3-3/10X

Résultats

Voici les formules des sommes calculées et factorisées à la fois à la main et par ordinateur pour les 8 premières valeurs de k :













On voit des relations entre toutes ces formules à deux indices d'écart (P2 est contenue dans P4, P6 et P8, P3 est contenue dans P5 et P7) mais je ne sais rien en faire pour l'instant. A poursuivre donc. on peut disposer les coefficients des polynômes successifs dans un tableau pour y voir plus clair, en y ajoutant quelques lignes supplémentaires pour donner plus de poids à nos conjectures :

kCoeff.012345678910111213
Les lignes 9 à 12 contiennent des coefficients calculés par le programme de cette page et n'ont pas été obtenus à la main...

On peut faire les constats suivants :

  • le coefficient de degré le plus élevé est toujours égal à 1 :
  • celui de plus bas degré est toujours nul, i.e. pas de terme constant ce qui est assez logique car :
  • la diagonale reliant les coefficients de degré k est une série arithmétique de raison 1/2 :
  • la diagonale située en-dessous est également une série avec pour formule :
  • quand la colonne 1 contient 0 (ce qui semble être le cas une fois sur deux) la diagonale issue de ce 0 ne contient que des 0 :
  • les autres diagonales sont également des séries avec pour formule :
  • La somme des coefficients de la ligne k est égale à k+1, ce qui permet de trouver les par simple soustraction.
Tout ceci semble indiquer qu'il y a une sorte de structure d'une ligne à l'autre. Je dis bien semble car je peux montrer comment calculer une ligne à partir de celle qui la précède mais je ne l'ai pas encore démontré rigoureusement, seulement vérifié jusqu'à un certain indice, c'est loin d'être la même chose.demonstration needed

� propos de la dernière remarque concernant les coefficients de la première colonne non-nulle je peux effectivement me débrouiller pour les calculer car la somme des autres coefficients est égale à k+1, mais il est important de remarquer que multipliés par le coefficient 1/k+1 ils sont égaux aux nombres de Bernouilli ! On vient donc de retrouver une méthode pour calculer ces nombres analogue à celle du triangle de Pascal pour les coefficients du binôme !formula needed.

Ce qui me fait croire qu'une partie de la démonstration des formules précédentes et de la procédure utilisée pour les obtenir pourrait résider dans ces fameux nombres de Bernouilli mais ça reste un terrain à investiguer.

Usages

Ces formules sont des ingrédients que l'on peut utiliser dans une décomposition polynomiale des termes d'une série (citation neeeded :)).

Questions ouvertes

Qu'en est-il des cas avec k non entier ? A faire ! �a peut sembler intéressant de formuler la sommes des racines carrées des N premiers entiers non ?

Références

Un article complet sur les sommes
La fameuse formule de Faulhaber
L'article de Knuth sur Faulhaber
Un mathématicien/logisticien/généticien qui s'est intéressé à la formule de Faulhaber
Les nombres de Bernouilli

Accueil1 2 3