Discussion :

[nakor]

Bonjour,
je cherche une manière de résoudre numériquement erf(x) qui est la fonction erreur, défini par

(http://fr.wikipedia.org/wiki/Fonction_d%27erreur)

Je ne veux pas utiliser directement les solutions existantes de la fonctions, mais je suppose qu'un algorithme permet de résoudre cette fonction. Peut être en utilisant son développement en série entière ? une idée ?

Merci

[FR119492]

Salut!
Voici un code Fortran qui fait exactement ce dont tu as besoin. Il ne te reste plus qu'à le décortiquer et à le traduire dans le langage de ton choix.

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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
      FUNCTION DERF(Y)
C                                  SPECIFICATIONS FOR ARGUMENTS
      DOUBLE PRECISION   DERF,Y
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      DIMENSION          P(5),Q(4),P1(9),Q1(8),P2(6),Q2(5)
      DOUBLE PRECISION   P,Q,P1,Q1,P2,Q2,XMIN,XLARGE,SQRPI,X,
     *                   RES,XSQ,XNUM,XDEN,XI
      INTEGER            ISW,I
C                                  COEFFICIENTS FOR 0.0 .LE. Y .LT.
C                                  .477
      DATA               P(1)/113.8641541510502D0/,
     *                   P(2)/377.4852376853020D0/,
     *                   P(3)/3209.377589138469D0/,
     *                   P(4)/.1857777061846032D0/,
     *                   P(5)/3.161123743870566D0/
      DATA               Q(1)/244.0246379344442D0/,
     *                   Q(2)/1282.616526077372D0/,
     *                   Q(3)/2844.236833439171D0/,
     *                   Q(4)/23.60129095234412D0/
C                                  COEFFICIENTS FOR .477 .LE. Y
C                                  .LE. 4.0
      DATA               P1(1)/8.883149794388376D0/,
     *                   P1(2)/66.11919063714163D0/,
     *                   P1(3)/298.6351381974001D0/,
     *                   P1(4)/881.9522212417691D0/,
     *                   P1(5)/1712.047612634071D0/,
     *                   P1(6)/2051.078377826071D0/,
     *                   P1(7)/1230.339354797997D0/,
     *                   P1(8)/2.153115354744038D-8/,
     *                   P1(9)/.5641884969886701D0/
      DATA               Q1(1)/117.6939508913125D0/,
     *                   Q1(2)/537.1811018620099D0/,
     *                   Q1(3)/1621.389574566690D0/,
     *                   Q1(4)/3290.799235733460D0/,
     *                   Q1(5)/4362.619090143247D0/,
     *                   Q1(6)/3439.367674143722D0/,
     *                   Q1(7)/1230.339354803749D0/,
     *                   Q1(8)/15.74492611070983D0/
C                                  COEFFICIENTS FOR 4.0 .LT. Y
      DATA               P2(1)/-3.603448999498044D-01/,
     *                   P2(2)/-1.257817261112292D-01/,
     *                   P2(3)/-1.608378514874228D-02/,
     *                   P2(4)/-6.587491615298378D-04/,
     *                   P2(5)/-1.631538713730210D-02/,
     *                   P2(6)/-3.053266349612323D-01/
      DATA               Q2(1)/1.872952849923460D0/,
     *                   Q2(2)/5.279051029514284D-01/,
     *                   Q2(3)/6.051834131244132D-02/,
     *                   Q2(4)/2.335204976268692D-03/,
     *                   Q2(5)/2.568520192289822D0/
C                                  CONSTANTS
      DATA               XMIN/1.0D-10/,XLARGE/6.375D0/
      DATA               SQRPI/.5641895835477563D0/
C                                  FIRST EXECUTABLE STATEMENT
      X = Y
      ISW = 1
      IF (X.GE.0.0D0) GO TO 5
      ISW = -1
      X = -X
    5 IF (X.LT..477D0) GO TO 10
      IF (X.LE.4.0D0) GO TO 25
      IF (X.LT.XLARGE) GO TO 35
      RES = 1.0D0
      GO TO 50
C                                  ABS(Y) .LT. .477, EVALUATE
C                                  APPROXIMATION FOR DERF
   10 IF (X.LT.XMIN) GO TO 20
      XSQ = X*X
      XNUM = P(4)*XSQ+P(5)
      XDEN = XSQ+Q(4)
      DO 15 I=1,3
         XNUM = XNUM*XSQ+P(I)
         XDEN = XDEN*XSQ+Q(I)
   15 CONTINUE
      RES = X*XNUM/XDEN
      GO TO 50
   20 RES = X*P(3)/Q(3)
      GO TO 50
C                                  .477 .LE. ABS(Y) .LE. 4.0
C                                  EVALUATE APPROXIMATION FOR DERF
   25 XSQ = X*X
      XNUM = P1(8)*X+P1(9)
      XDEN = X+Q1(8)
      DO 30 I=1,7
         XNUM = XNUM*X+P1(I)
         XDEN = XDEN*X+Q1(I)
   30 CONTINUE
      RES = XNUM/XDEN
      GO TO 45
C                                  4.0 .LT. ABS(Y), EVALUATE
C                                  MINIMAX APPROXIMATION FOR DERF
   35 XSQ = X*X
      XI = 1.0D0/XSQ
      XNUM = P2(5)*XI+P2(6)
      XDEN = XI+Q2(5)
      DO 40 I=1,4
         XNUM = XNUM*XI+P2(I)
         XDEN = XDEN*XI+Q2(I)
   40 CONTINUE
      RES = (SQRPI+XI*XNUM/XDEN)/X
   45 RES = RES*DEXP(-XSQ)
      RES = 1.0D0-RES
   50 IF (ISW.EQ.-1) RES = -RES
      DERF = RES
      RETURN
      END

Jean-Marc Blanc