Répondre

yoshi · 08-01-2013 10:30:33

Bonjour,

Toujours risqué de passer après (et encore pire, avant !)...
Bon, au départ je me suis mépris sur Scilab, puis j'ai découvert l'avoir déjà installé...
Alors j'ai testé le petit bout de programme proposé.
1. Il fonctionne tel quel (à la modification de fin de boucle près), la preuve :
Pour p = 10

2. Le calcul est-il exact ?
Résultat Python :

IDLE 2.6.4
>>> ================================ RESTART ================================
>>>
S = 0.995068212238
>>>

Donc oui.
Le calcul direct de l'intégrale est-il exact (pourquoi non ?) ?
Résultat avec le logiciel de calcul formel libre et gratuit WxMaxima :
0.995054031001174

Oui aussi.

Alors, que se passe-t-il ?

Pour le plaisir, 2 versions Python, l'une "bête et méchante" (!) :

# usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division

from math import sqrt

p,n=10,100

x0,x4=1,10

h=(x4-x0)/n

S=0

def f(a,p):

    y=1/(sqrt(1+(a-1/a)**2*p))

    return y

for nb in xrange(0,n):

    x=x0+nb*h

    S+=(f(x,p)+f(x+h,p))*h/2

print "S =",S

l'autre "factorisée" et demandant le minimum d'opérations au programme. :

# usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division

from math import sqrt

p,n=10,100

x0,x4=1,10

h,S=(x4-x0)/n,0

def f(a,p):

    return 1/(sqrt(1+(a-1/a)**2*p))

for nb in xrange(1,n):

    S+=f(x0+h*nb,p)

S=((f(x0,p)+f(x4,p))/2+S)*h

print "S =",S

Explications sur "factorisée" :
Si a0 = f(x0), a1 =f(x0+h), a2 = f(x0+2h), a3 = f(x0+3h)...an=f(x0+n*h)
le calcul de la somme s'écrit :
[tex]S=(a_0+a_1)\times \frac h 2 + (a_1+a_2)\times \frac h 2 + (a_2+a_3)\times \frac h 2 +... +(a_{n-1}+a_n)\times\frac h 2[/tex]

Soit [tex]S = [(a_0+a_1)+(a_1+a_2)+(a_2+a_3)+...+(a_{n-2}+a_{n-1})]\times \frac h 2[/tex]

Soit encore [tex]S=(a_0+2a_1+2a_2+2a_3+...+2a_{n-1}+a_n)\times \frac h 2[/tex]

Et enfin :

[tex]S= \left(\frac{a_0+a_n}{2}+a_1+a_2+...+a_{n-1}\right)\times h[/tex]
La boucle se contentera donc de calculer [tex]\sum_{i=1}^{n-1} a_i[/tex]
Ça allait peut-être sans dire, mais ça allait encore bien mieux en l'explicitant...

C'est ainsi qu'on demande le moins de calculs au processeur.

@+

totomm · 07-01-2013 15:28:37

Bonjour,

En SciLab le for est défini : for variable = début : pas : fin, …, end
Il faut donc sans doute écrire :
for x=x0:h:x4-h; plutôt que for x=x0:h:x4; pour n'avoir que n évaluation des trapèzes
(le premier de x0 à x0+h et le nième de (x0+(n-1)h = x4-h) à (x0+nh) = x4

En programmation de la méthode des trapèzes mieux vaut évaluer la demi-somme de la fonction pour x_début et x_fin plus la somme pour les x entre x_début et x_fin
Il faut aussi se méfier des pas exprimés en float ou double : exprimer les pas par des entiers, en ayant une variable x qui s'incrémente en float ou double...

pour comparaison voici mon estimation de votre intégrale telle que reprise en post #1 par yoshi
pour différentes valeurs de p et n=100
p = 0.5 : S = 3.12326
p = 1 : S = 2.44392
p = 2 : S = 1.88608
p = 5 : S = 1.31754

Cordialement

factions31 · 06-01-2013 22:42:05

bonjour, j'ai besoin d'un peu d'aide si possible

Algorithme pour le calcul de l’erreur avec la méthode des trapèzes

clear
//calcul des valeurs approximatives de l'intégrale
n=input('donner la valeur de segmentations n ');

x0=1
x4=10
h=(x4-x0)/n;
x=x0:h:x4;
p=input('donner la valeur de p')

function y=f(x)
y=1/(sqrt(1+p*(x-(1/x))^2));
endfunction
S=0
for x=x0:h:x4;
S=S+h*(f(x)+f(x+h))/2;
end
S

//calcul de la valeur exacte de l’intégrale
i=integrate('1/(sqrt(1+p*(x-(1/x))^2))','x',x0,x4)

//calcul de l’erreur
i-S

cet algorithme ne marche pas et je ne sais pas pourquoi quelqu'un pourrait me dire où est mon erreur?

merci d'avance :D

Forum de mathématiques - Bibm@th.net

Répondre

Résumé de la discussion (messages les plus récents en premier)

Pied de page des forums