Détermination numérique temps de relaxation

[invite9d8ab2d1]

Bonjour,

J'ai une équation différentielle de la forme d²y/dy² + 1/Qdy/dx + y = 0 et je doit trouver numériquement le temps de relaxation. C'est à dire le temps t_0 à partir duquel pour tous les t > à t_0, abs(y(t)) < 1/100.
J'ai donc résolu l'équation et j'ai les solutions et la courbe de tracé (en Python).
Là où je sèche c'est sûr la façon de déterminer t_0. J'avais pensé à une boucle while toute simple mais lorsque Q > 1/2, on est en régime pseudo périodique donc la fonction devient inférieur à 1/100 plusieurs fois comme on voit ici :
Figure_1.png
Mais la courbe peut aussi ressembler à ça :
Figure_1.png

Après peut être que je pourrai trouver la courbe qui suit les maximums et déduire t_0 de cela mais je ne crois pas avoir vu ça en cours.
Voilà après plusieurs jours de recherche sur internet je bloque et donc si vous avez une piste je suis preneur !

Merci
-----

[invite6c250b59]

Une façon simple est de prendre un critère sur une fenêtre glissante, genre max(t:t+delta) - min(t:t+delta)

[invite23cdddab]

Plus simple : si tu as les solutions numériques, tu as un tableau de valeurs de y(t).

Donc tu parcours le tableau, et si |y(t_n)| >= 1/100, alors tu actualise la valeur de T_relaxation à la valeur t_n.

Une fois que tu as parcouru tout le tableau, T_relaxation est la valeur que tu recherches (puisque c'est le dernier moment ou |y(t)|> 1/100 )

[pm42]

Puisqu'on a le tableau pour tracer la courbe, on peut en effet chercher en partant de la fin le 1er élément qui est plus grand que 1/100.

Cela donne :

next(i for i, e in reversed(list(enumerate(tableau))) if abs(e) > 0.01)+1

en supposant que tableau contient les valeurs. On peut le coder d'une façon plus longue mais c'est moins rigolo.

[A voir en vidéo sur Futura]

[invite6c250b59]

Roo c'est joli! Allez pinaillage quand même

next(i for i, e in reversed(list(enumerate(tableau, 1))) if abs(e) > 0.01)

[pm42]

Roo c'est joli! Allez pinaillage quand même

Joli pinaillage aussi

[invite9d8ab2d1]

Merci beaucoup à tous, j'ai essayé la méthode avec le tableau de valeur et ça fonctionne bien c'est parfait !
Il faut juste bien choisir la précision pour créer le tableau

[invite6c250b59]

Petite exploration de l'idée d'utiliser un itérateur.

import math         
myfunc = lambda t : math.pi/t**1/10
        
class convergence:
    # allows computing for which t your function goes below your required precision 
    # beware too large precision and non converging function

    def __init__(self, func, precision, verbose):
        tf, precision = 1, 1/10**precision
        while func(tf) > precision:
            ti, tf = tf, tf*2
        self.ti = ti
        self.tf = tf
        self.p = precision
        self.f = func
        self.v = verbose
        
    def __iter__(self):
        return self
            
    def __next__(self):
        tm = (self.tf + self.ti) / 2
        if self.f(tm)>self.p:
            self.ti = tm
        else:
            self.tf = tm
        
        if self.v:
            print(self.ti, self.tf)
        
        return self.tf-self.ti>self.p

_ = convergence(myfunc, 5, True)
while next(_): pass

print("")
print("convergence pour t = ", _.tf)
print("avec f(t) = ", _.f(_.tf))

[invite9dc7b526]

avec R, en supposant que x est le vecteur contenant les valeurs de la fonction

(length(x)+1)-(which(cumsum(x[length(x):1]>0.01)==1)-1)

[invite6c250b59]

Avec matlab, en supposant que x contient la solution:

x

[invite6c250b59]

Petite exploration de l'idée d'utiliser un itérateur.

Pour revenir là-dessus, au final l'intérêt n'est pas si évident (ou du moins pas de la manière que j'ai essayée).

def myf(func, precision, verbose):
        tf, precision = 1, 1/10**precision
        while func(tf) > precision:
            ti, tf = tf, tf*2
        
        while tf-ti>precision:
            tm = (tf + ti) / 2
            if func(tm)>precision:
                ti = tm
            else:
                tf = tm
            if verbose:
                print(ti, tf)        
        return tf, func(tf)

a, b = myf(myfunc, 8, True)