Intégration numérique d'une fonction

**Helber** · 17/05/2010, 13h41

Bonjour,

Je suis presque nouveau sur Delphi mais j'y consacre en ce moment pas mal de temps.
J'aurais besoin pour mon projet d'avoir une fonction qui me permettrait d'intégrer (au sens calcul de l'intégrale entre 2 bornes) une autre fonction.
Dans l'idéal cette fonction "intégration" serait définie dans un fichier à part. Ensuite si dans un autre fichier j'ai besoin d'intégrer une fonction, j'aurais juste à appeler ce fichier où se trouve cette fonction intégration et à écrire à l'endroit voulu dans ma procédure " mavariable:=integration(mafonctionàintégrer,borne_inf,borne_sup);". J'ai lutté ces derniers jours pour y arriver ; sans succès.
Ma recherche sur le forum n'a pas donné plus de résultat...
est-ce que par hasard quelqu'un a déjà implémenté cet outil? Ou a une idée assez précise de comment on peut faire?
Merci par avance,
Helber

**Charly910** · 17/05/2010, 14h36

Bonjour,

essaye ceci (à perfectionner) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
unit Main;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;
type
  TOperation = function(const X: Extended): Extended;
 
type
  TF_Main = class(TForm)
    Btn_Integration: TButton;
    procedure Btn_IntegrationClick(Sender: TObject);
  private
    { Déclarations privées }
  public
    { Déclarations publiques }
  end;
 
var
  F_Main: TF_Main;
 
implementation
 
{$R *.dfm}
 
{ ========================================================================= }
function Integrale(const X1,X2 : Extended; const Fonction: TOperation): Extended;
Var
   N, i : Integer ;
   Somme, XDeb, XFin, Pas : Extended ;
   YDeb, YFin : Extended ;
begin
  N := 100 ;
  Somme := 0 ;
  Pas := (X2 - X1)/N ;
  For i := 1  to N Do
     Begin
        XDeb := X1 + (i-1) * Pas ;
        XFin := Xdeb + Pas ;
        YDeb :=  Fonction(XDeb) ;
        YFin :=  Fonction(XFin) ;
        Somme := Somme + (YDeb + YFin) / 2 * Pas ;
     End ;
  result := Somme ;
end;
{ ========================================================================= }
function Carre(const X: Extended): Extended ;
begin
  result := X * X;
end;
{ ========================================================================= }
function Cube(const X: Extended): Extended ;
begin
  result := X * X * X ;
end;
{ ========================================================================= }
procedure Test;
var X,Y,A,M : Extended ;
begin
  X := 1.00 ;
  Y := 2.00 ;
  A := Integrale(X, Y, Carre);
  M := Integrale(X, Y, Cube);
//  ShowMessage(Format('X= %*.*f, Y= %*.*f // Intégrale Carré= %*.*f // Intégrale Cube= %*.*f',[X, Y, A, M]));
  ShowMessage ('Intégration entre :'+Chr(10)+Chr(13)+
               '    X1 = '+FloatToStr(X)+Chr(10)+Chr(13)+
               '    X2 = '+FloatToStr(Y)+Chr(10)+Chr(13)+
               '       Fonction y=x2 :  '+FloatToStr(A)+Chr(10)+Chr(13)+
               '       Fonction y=x3 :  '+FloatToStr(M));
 
 
  end;
{ ========================================================================= }
procedure TF_Main.Btn_IntegrationClick(Sender: TObject);
begin
   Test ;
end;
 
end.

tu mets le bouton "Btn-Integration" qui appelle la fonction Test dans ta form

A+

Charly

**Helber** · 17/05/2010, 15h07

Merci pour cette première réponse. J'ai essayé de faire quelque chose avec mais ça n'a pas l'air de marcher

.
J'ai écrit la fonction qui intègre dans une Unit "Integration" (fichier .pas) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
unit Integration;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;
type
  TOperation = function(const X: Extended): Extended;
 
implementation
 
   function Integrale(const X1,X2 : Extended; const Fonction: TOperation): Extended;
Var
   N, i : Integer ;
   Somme, XDeb, XFin, Pas : Extended ;
   YDeb, YFin : Extended ;
begin
  N := 100 ;
  Somme := 0 ;
  Pas := (X2 - X1)/N ;
  For i := 1  to N Do
     Begin
        XDeb := X1 + (i-1) * Pas ;
        XFin := Xdeb + Pas ;
        YDeb :=  Fonction(XDeb) ;
        YFin :=  Fonction(XFin) ;
        Somme := Somme + (YDeb + YFin) / 2 * Pas ;
     End ;
  result := Somme ;
end;
 
end.

Ensuite, dans une autre unit déjà existante (aussi fichier .pas, dans le même dossier où j'ai mis "Integration"), j'ai appelé "Integration" dans "uses" et ai utilisé la fonction integrale définie dans Integration comme suit :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
unit DGL1D_Interface_ThermionicEmission;

interface

uses Dialogs,SysUtils,DGL1D_Interface,grdpnt,Direction_Parent, Math,Integration;

//type TFunctionParam = function(x:extended):extended;

type
// Semiconductor-Semiconductor Thermionic-Emission Interface
TDGL1D_Interface_ThermionicEmission =  class(TDGL1D_Interface)

XXX

  // DGLsystem to be solved
  procedure DGLSys; override;                      //DC mode DGL

end;


implementation


uses HeteroJunctionMain,type_def,ParameterIO,constants,Direction_Interface,
     Direction_Interface_TE,ModelParametersTE;


    function test(B: extended): extended;
begin
  result := B*5;
end;

//steady-state DC DGLSys, Thermionic Emission Interface
procedure TDGL1D_Interface_ThermionicEmission.DGLSys;


var

  inte: extended;


begin

inte:=Integrale(1.0,2.0,test);

end;

Bon bien sûr il y a pas mal d'autre ligne de code dans la version originale mais pour alléger j'ai supprimé ce qui ai extérieur au problème.
Donc où est l'erreur maintenant?! Faut-il déclarer "Integration" ailleurs?
Merci encore!
Helber

**Stevens.Georges** · 17/05/2010, 15h15

Bonjour,

A ma connaissance il n’est pas possible de disposer d’une fonction (Delphi) simple et générale retournant la valeur d’une intégrale entre les limites inférieure et supérieure. En effet, cette fonction (Delphi) devrait intégrer la fonction (mathématique) dont on cherche l’intégrale, or il y en a évidemment une kyrielle.

J’ai eu personnellement quelques problèmes de ce type à résoudre et y suis parvenu en suivant le processus utilisé par les calculettes scientifiques par exemple HP 33 S dont tu trouveras ci-dessous la description sommaire.

Li : limite inférieure.
Ls : limite supérieure.
f(x) : la fonction qu’on cherche à intégrer.

1°/ Diviser la plage Ls – Li en 5 tranches.
2°/ Pour chaque tranche calculer la valeur moyenne Xm1, Xm2, Xm3, Xm4, Xm5.
3°/ Calculer f(x) pour chaque valeur de Xm.
4°/ Totaliser T5 = f(Xm1) + f(Xm2) + f(Xm3) + f(Xm4) + f(Xm5) .

Recommencer exactement selon le même processus en doublant le nombre de tranches soit 10 et comparer T5 et T10.

Et cetera.

Arrêter lorsque la différence entre TN et TN-1 devient négligeable.

Ca n’est évidemment pas rigoureusement exact au sens mathématique du terme mais en choisissant intelligemment la valeur de la différence TN – TN-1 ça marche parfaitement.

**Helber** · 17/05/2010, 15h26

Bonjour Stevens.Georges,

il me semble que tu confonds primitive et intégrale : si en effet il y a une kyrielle de primitive pour une fonction donnée, il n'y a qu'une seule intégrale de cette fonction entre 2 bornes (la constante additive s'annule si c'est à ça que tu penses).
Mon problème n'est pas au niveau mathématique mais bien au niveau de "comment inclure ça dans Delphi". Après il existe pleins d'algorithmes plus ou moins compliqués et précis donnant le bon résultat (dont celui que tu donnes!).
Si tu pouvais m'indiquer comment tu as implémenté cette fonction dans ton code (je veux dire, la structure générale) ça serait super!
Merci
Helber

**Charly910** · 17/05/2010, 15h39

@ Stevens.Georges

c'est vrai que l'on ne peut calculer qu'une approximation comme dans mon exemple par la méthode des trapèzes (j'ai pris 100 intervalles, pour plus de précision on peut augmenter n) - on peut aussi prendre la méthode de Simpson qui est un peu plus précise.

@ Helber

je pense que cela ne fonctionne pas car tu n'as pas déclaré la fonction Integrale dans le module Integration

Voici mon exemple coupé en 2 modules

L'unité Intégration :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
unit Integration;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;
type
  TOperation = function(const X: Extended): Extended;
 
  function Integrale(const X1,X2 : Extended; const Fonction: TOperation): Extended;
 
 
implementation
 
   function Integrale(const X1,X2 : Extended; const Fonction: TOperation): Extended;
Var
   N, i : Integer ;
   Somme, XDeb, XFin, Pas : Extended ;
   YDeb, YFin : Extended ;
begin
  N := 1000 ;
  Somme := 0 ;
  Pas := (X2 - X1)/N ;
  For i := 1  to N Do
     Begin
        XDeb := X1 + (i-1) * Pas ;
        XFin := Xdeb + Pas ;
        YDeb :=  Fonction(XDeb) ;
        YFin :=  Fonction(XFin) ;
        Somme := Somme + (YDeb + YFin) / 2 * Pas ;
     End ;
  result := Somme ;
end;
 
end.

et le module principal :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
unit Main;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, Integration ;
type
  TOperation = function(const X: Extended): Extended;
 
type
  TF_Main = class(TForm)
    Btn_Integration: TButton;
    procedure Test;
    procedure Btn_IntegrationClick(Sender: TObject);
  private
    { Déclarations privées }
  public
    { Déclarations publiques }
  end;
 
var
  F_Main: TF_Main;
 
implementation
 
{$R *.dfm}
 
{ ========================================================================= }
function Carre(const X: Extended): Extended ;
begin
  result := X * X;
end;
{ ========================================================================= }
function Cube(const X: Extended): Extended ;
begin
  result := X * X * X ;
end;
{ ========================================================================= }
procedure TF_Main.Test;
var X,Y,A,M : Extended ;
begin
  X := 1.00 ;
  Y := 2.00 ;
  A := Integrale(X, Y, Carre);
  M := Integrale(X, Y, Cube);
  ShowMessage ('Intégration entre :'+Chr(10)+Chr(13)+
               '    X1 = '+FloatToStr(X)+Chr(10)+Chr(13)+
               '    X2 = '+FloatToStr(Y)+Chr(10)+Chr(13)+
               '       Fonction y=x2 :  '+FloatToStr(A)+Chr(10)+Chr(13)+
               '       Fonction y=x3 :  '+FloatToStr(M));
  end;
{ ========================================================================= }
procedure TF_Main.Btn_IntegrationClick(Sender: TObject);
begin
   Test ;
end;
 
end.

Il faut aussi déclarer Integration dans le Uses du Main

Charly

**Helber** · 17/05/2010, 15h56

Ahah! ça va de mieux en mieux! J'ai fait ta correction suggérée , puis enlever les "const" dans les déclarations des fonctions (à quoi servent-ils?!) car mes fonctions à intégrer n'en avait pas. Et maintenant ça marche... à un petit détail près : seule la fonction définie en dehors de "procedure" veut bien se faire intégrer : ça marche avec "test", mais pas avec "theta"...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
unit DGL1D_Interface_ThermionicEmission;

interface

uses Dialogs,SysUtils,DGL1D_Interface,grdpnt,Direction_Parent, Math,Integration;


type
// Semiconductor-Semiconductor Thermionic-Emission Interface
TDGL1D_Interface_ThermionicEmission =  class(TDGL1D_Interface)

  vari:extended;

  // DGLsystem to be solved
  procedure DGLSys; override;                      //DC mode DGL


end;


implementation


uses HeteroJunctionMain,type_def,ParameterIO,constants,Direction_Interface,
     Direction_Interface_TE,ModelParametersTE;

    function test(B: extended): extended;
begin
  result := B*5;
end;

//steady-state DC DGLSys, Thermionic Emission Interface

procedure TDGL1D_Interface_ThermionicEmission.DGLSys;


var

  inte: extended;

  //Heavyside function
  function theta(xx:Extended):Extended;
  begin
  result := 0;
  if xx<0 then result :=1;
  end;

begin
//initialisation
for i:=1 to anzDGL do for j:=1 to anzDgl do begin
  gp.DGL_0[i,j]:=0;
  gp.DGL_0Im[i,j]:=0;
  gp.DGL_F[i,j]:=0;
//???  gp.fmc[i,j]:=0;
  gp.DGL_B[i,j]:=0;
  end;

...

inte:=Integrale(1.0,2.0,theta);

end;

Qu'est-ce qui cloche maintenant?
Merci tout de même déjà pour cette avancée!

**Charly910** · 17/05/2010, 16h34

Je ne sais pas trop pourquoi Theta ne fonctionne pas

c'est surement une question de visibilité de la fonction integrale (un pro du forum pourrait surement te l'expliquer)

Est ce que cela te gêne de mettre theta hors de la procédure ?

Sinon pour les Const, c'est une question de déclaration :

si un paramètre est en entrée seulement on met const, sinon on met var, mais ce n'est pas obligatoire, tu peux enlever les const.

A+

Charly

**Stevens.Georges** · 17/05/2010, 16h58

Je n'ai pas confondu intégrale et primitive. J'ai seulement voulu dire qu'une function Delphi universelle devrait connaitre toutes (kyrielles) les fonctions à intégrer.

Voici comment je procède, sauf erreur de ma part car j'ai re-écrit le code le mien réel étant trop complexe pour être intelligible.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
implementation
 
{$R *.DFM}
 
var N:integer; //Nombre de tranches
    TN:single;//f(Xm1)+f(Xm2)+f(Xm3)...+f(XmN)
    TNm1:single;//f(Xm1)+f(Xm2)+f(Xm3)...+f(XmN-1)
    Li,Ls:single;//Limites inférieure supérieure
    TR:single;//Tranche
    DifLim:single;//Différence limite.
    Dif,Xm:single;
 
function F(X:single):single;
begin
  //Calcule la valeur de la fonction f dont
  //on fait l'intégration pour x.
end;
 
procedure Integration;
var I:integer;
begin
  TNm1:=0;
  repeat
    Tr:=(Ls-Li)/N;TN:=0;
    for I:=1 to N do
      begin
        Xm:=(N-0.5)*Tr;//Déterminer Xm en milieu de tranche
        TN:=TN+F(Xm);
      end;
    Dif:=TN-TNm1;
    if abs(Dif)>DifLim then
      begin //Si la différence > différence limite
            //augmenter N et retenir TN comme TNm1
        Inc(N);TNm1:=TN;
      end;
  until abs(Dif)<=DifLim;
end;

**Stevens.Georges** · 18/05/2010, 08h59

Bonjour,

Je me suis aperçu d'une omission.

La ligne ci-dessous :

TN:=TN+F(Xm);

doit être corrigée comme suit:

TN:=TN+F(Xm)*Tr;

A+

**Helber** · 18/05/2010, 09h10

Bonjour,

Merci pour votre aide. Il me reste juste ce problème de non visibilité de ma fonction à intégrer lorsqu'elle est défini dans ma procédure. Est-ce qu'un pro du forum passe par là pour éventuellement m'expliquer ce problème? (qui pourra me gêner assez rapidement je pense...)
Helber

**Franck SORIANO** · 19/05/2010, 20h57

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
function Integrale(const X1,X2 : Extended; const Fonction: TOperation): Extended;
Var
   N, i : Integer ;
   Somme, XDeb, XFin, Pas : Extended ;
   YDeb, YFin : Extended ;
begin
  N := 1000 ;
  Somme := 0 ;
  Pas := (X2 - X1)/N ;
  For i := 1  to N Do
     Begin
        XDeb := X1 + (i-1) * Pas ;
        XFin := Xdeb + Pas ;
        YDeb :=  Fonction(XDeb) ;
        YFin :=  Fonction(XFin) ;
        Somme := Somme + (YDeb + YFin) / 2 * Pas ;
     End ;
  result := Somme ;
end;

Désolé, mais ça c'est une très mauvaise façon de calculer l'intégrale... Ca peut fonctionner avec certaines fonctions, mais on risque aussi d'obtenir des résultats très aléatoires en fonction de la fonction à intégrer et de l'amplitude de l'intervale d'intégration. Autrement dit, ça ne marche pas dans un cas général.

En effet, pour faire l'intégration numérique, on utilise ici la méthode des trapèzes.
Cependant, celà revient à échantillonner la fonction que l'on veut intégrer.
On doit alors faire attention à vérifier le théorème de Nyquist-Shannon :
En un mot, il faut que la fréquence d'échantillonnage soit au moins égale au double de la fréquence maximale du signal échantillonné.
Sinon, on risque de faire apparaître un phénomène de repliement de spectre, et le résultat de l'échantillonnage ne reflétera pas la fonction échantillonnée.

Par exemple, avec cette implémentation :
Integrale(0, 2*pi, cos) = -3.09 E-20.
Normalement, le résultat exact devrait être 0. On a obtenu un résultat proche de 0. On peut éventuellement accepter l'erreur.

Par contre :
Integrale(0, pi*2000, cos) = 6283
Cette fois c'est du grand n'importe quoi... cos(x) étant périodique, on aurait du avoir 1000*Integrale(0, 2*pi, cos), soit encore 0, au lieu de ça on est tombé sur 2*pi*1000 !

**Helber** · 20/05/2010, 11h34

Oui oui tout à fait d'accord : l'algorithme est vraiment simpliste mais je veux avoir la structure générale avant de regarder l'aspect mathématique des choses!

**Stevens.Georges** · 21/05/2010, 15h47

C'est évident qu'il y a des précautions à prendre avant de lancer l'intégration.

Quant à dire que c'est une très mauvaise méthode, je veux bien mais c'est celle qu'utilise ma HP 33S. Y a-t-il d'ailleurs une meilleure méthode sauf bien sur à connaitre la fonction à intégrer et sa primitive.

**Charly910** · 22/05/2010, 10h25

Bonjour,

il n'y a pas de meilleure méthode d'intégration numérique (sauf Simpson qui améliore un peu le résultat car il prend une parabole passant par 3 points successifs au lieu d'une droite entre 2 points pour la méthode des trapèzes)

la précision obtenue est fonction du Pas [ (xb - xa ) / n ] ; l'intégrale est juste quand le pas tend vers 0 ( c'est à dire quand n tend vers l'infini)

c'est pourquoi le 2ième exemple de Franck ne donne pas de bon résultat : pas beaucoup trop grand

On ne peut vérifier la validité du résultat qu'en faisant augmenter n et en testant la convergence du résultat comme l'a proposé Stevens.Georges

A+

Charly

**Franck SORIANO** · 23/05/2010, 18h30

Envoyé par Stevens.Georges

Quant à dire que c'est une très mauvaise méthode, je veux bien mais c'est celle qu'utilise ma HP 33S.

En fait je parlais de la solution de Charly910 qui utilise un pas fixe, calculé à partir de l'amplitude de l'intervale d'intégration.
Comme l'a indiquée Charly910 cette méthode fonctionne lorsque le pas tend vers 0. Le problème c'est qu'il faut avoir un pas adapté à la fonction.

Intégration numérique d'une fonction

Langage Delphi

Discussions similaires

Partager

Partager