Forum de mathématiques - Bibm@th.net

factions31

Invité

scilab erreur

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

totomm

Membre

Inscription : 25-08-2011

Messages : 1 093

Re : scilab erreur

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

Dernière modification par totomm (07-01-2013 15:35:25)

yoshi

Modo Ferox

Inscription : 20-11-2005

Messages : 17 385

Re : scilab erreur

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.

@+