Calcul d’intégrale fractionnaire Python

cristiano_r
cristiano_r
dans Logiciels pour les mathématiques
Bonjour à tous
Voici mon problème, j'aimerais calculer des intégrales du type $$
\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydz\vert x-y\vert ^{H-1} \vert y-z\vert ^{H-1}\vert z-x\vert ^{H-1}
$$ à l'aide de Python, je connais la valeur exacte de celle-ci qui s'exprime avec la fonction Beta (l'objectif étant de pouvoir changer les bornes avec des $t_1,t_2,t_3$ différents où il n'y a pas forcément de formule explicite).

Cependant en utilisant la fonctions "scipy.integrate.tplquad", il y a déjà le problème d'avoir "0.0 cannot be raised to a négative power", j'ai donc rajouté des epsilon très petits dans la fonction elle même, mais Python me dit que l’intégrale est probablement divergente ou converge lentement.

J'ai donc codé bêtement les sommes de Riemann sans aucune optimisation, là je peux avoir un résultat mais il est extrêmement sensible à la valeur d'epsilon dans la fonction et dans les bornes car commencer à partir de 0 renvoie aussi des erreurs.
Merci d'avance pour votre aide.

Réponses

  • Rescassol
    Rescassol
    Bonjour,

    Quelque chose comme ça ? 
    from scipy.integrate import nquad

def beta():
    def f(x,y,z):
        return abs(((x-y)*(y-z)*(z-x)))**(H-1)

    def r1(x):
        return [0,1]

    def r2(x,y):
        return [0,1]

    return nquad(f,[r2,r1,[0,1]])[0]

H=2.0
B=beta()
print(B)

    Cordialement,

    Rescassol
  • cristiano_r
    cristiano_r
    J’ai oublié de préciser que H est dans (1/2,1) les deux bornes exclus, merci pour le code j’avais la même idée cependant c’est ici le paramètre H qui est problématique
  • cristiano_r
    cristiano_r
    Je suis vraiment bloqué, je ne trouve pas de méthode alternative que des rajouter des epsilons qui faussent le résultat. Si quelqu’un à une idée du problème que Python rencontre avec ce type d’intégrale.
  • Rescassol
    Rescassol
    Bonjour

    Tu pourrais donner les valeurs attendues, ou alors la formule avec la fonction Beta ? Pour $H \in [\frac{1}{2}; 1]$.

    Cordialement,
    Rescassol
  • cristiano_r
    cristiano_r
    Bonsoir,

    Alors on trouve exactement $\frac{2}{H(3H-1)}Beta(H,H)$ soit avec H=0.75 environ 3.69. Ce que je n'arrive pas à retrouver avec Python sans passer par cette formule.
  • Rescassol
    Rescassol
    Bonsoir,

    Un peu de MonteCarlo: 
    def f(x,y,z):
    A=abs((x-y)*(y-z)*(z-x))
    return A**(H-1)

def betaMC():
    N=1000000
    cpt=0
    for k in range(N):
        T = np.random.rand(3)
        cpt += f(T[0],T[1],T[2])
    return cpt/N

H=0.75
B=betaMC()
print(B)

    Cordialement,

    Rescassol
  • cristiano_r
    cristiano_r
    Merci beaucoup Rescassol, ça fonctionne bien. Cependant j’ai toujours une question, sauriez-vous me dire pourquoi la fonction scripy.integrate n’arrive pas à calculer cette intégrale. Sinon encore merci.
  • Rescassol
    Rescassol
    Bonsoir,

    C'est parce que la fonction à intégrer n'est pas bornée sur le domaine d'intégration.
    Regarde ce qu'il se passe si au moins deux des trois $x,y,z$ sont égaux et $H<1$.

    Cordialement,

    Rescassol
Connectez-vous ou Inscrivez-vous pour répondre.

Bonjour!

Pour participer au forum, cliquer sur l'un des boutons :

Connexion S'inscrire

Catégories

In this Discussion

Qui est en ligne 17

AmadouBibixGaBuZoMeuLirone93NicoLeProfbiguine_equationzeste +10 Invités