Discussion :

[Razgriz]

Bonour à tous, voilà j'ai un problème avec RK4 :

Soit le problème de Cauchy

d u(t) = f(t, u(t) )
u(t0) = u0

Je veux trouver la fonction u(t) en intégrant avec Runge Kutta d'ordre 4 en intégrant sur l'intervalle [t0, t1]

Je crée donc une classe RungeKutta2D qui va faire ça. Seulement voilà, avec de bêtes exemples ça marche pas :

sur [0, 10] : d u(t) = t , u(0) = 0
==> je voudrais qu'on me retourne la fonction u(t) = t² / 2 car d(t² / 2) = t mis ça me donne pas ça

Je travaille en Java

Classe qui fait RungeKutta :

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
 
public class RungeKutta2D
{
    private double u0;
    private double t1;
    private double t0;
    private Function2D f;
 
    public RungeKutta2D(Function2D f, double t0, double t1, double u0)
    {
        this.f = f;
        this.t0 = t0;
        this.t1 = t1;
        this.u0 = u0;
    }        
 
    public Point2D.Double[] getRK4(int n)
    {
        Point2D.Double[] pts = new Point2D.Double[n];
        double ui = u0;
        pts[0] = new Point2D.Double(t0,u0);
        double h = (t1 - t0) / (n - 1);
 
        for (int i = 1; i < pts.length; i++)
        {
            double ti = t0 + i * h;
            double k1 = h * f.eval(ti, ui);
            double k2 = h * f.eval(ti + h/2 , ui + k1/2);
            double k3 = h * f.eval(ti + h/2 , ui+ k2/2);
            double k4 = h * f.eval(ti + h, ui + k3);
 
            ui += (1/6.) * (k1 + 2*k2 + 2*k3 + k4);
            pts[i] = new Point2D.Double(ti, ui);
        }
 
        return pts;
    }
}

Interface Function2D :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
public interface Function2D
{
            public double eval(double x, double y);
}

sur le bête test :

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
 
public class RKTest
{
 
    /**
     * @param args the command line arguments
     */
    public static void main(String[] args)
    {
        //si d u(t) = t alors f(t, u(t)) = t
        // => f : R² -> R : (x, y) ->  x
 
        //analytiquement, u(t) = x² / 2
 
        Function2D f = new Function2D()
        {
            public double eval(double x, double y)
            {
                return x;
            }
        };
 
        RungeKutta2D rk = new RungeKutta2D(f, 0, 10 , 0);
        Point2D.Double[] pts = rk.getRK4(5);
        for(Point2D.Double p : pts)
             System.out.println(p);

et ça m'imprime :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
Point2D.Double[0.0, 0.0]
Point2D.Double[2.5, 9.375]
Point2D.Double[5.0, 25.0]
Point2D.Double[7.5, 46.875]
Point2D.Double[10.0, 75.0]

or je devrais avoir (0,0) , (2.5 , 3.125) , (5 , 12.5) , (7.5 , 28.125) , (10, 50)
ou des trucs très proches tout du moins...

Qu'est-ce qui ne fonctionne pas???

[millie]

Regarde ma source :

http://www.developpez.net/forums/sho...d.php?t=236982

(edit, je croyais que c'était du C++ ta source (j'eusses regardé trop vite), mais ça doit quasiment être pareil )

[Razgriz]

Votre code me semble très correct. A vrai dire, tout aussi correct que le mien.

ma question est la suivante : pourquoi ça fonctionne pas avec l'exemple que je fais? Qu'est-ce que je fais qui faut pas faire (ou que je ne fais pas qu'il faut faire), où est la faute???

Merci d'avance, ça fait un moment que je cherche, et je trouve pas...

[millie]

Merci d'avance, ça fait un moment que je cherche, et je trouve pas...

Moi je sais !

Tu as n=5 et t0 = 0 et t1 = 10 soit un pas h de :

h = 10 / 4 ~ 2... Et donc, j'en déduis que le pas est beaucoup trop grand.

[Razgriz]

Ben non justement :
si t0 = 0 et t1 = 10 et que je fais des pas de (t0 - t1) / (n - 1) ça fait 10 / 4 = 2.5

j'aurais donc les 5 points (le 5 c'est le nombre de points) suivants :

0
2.5
5
7.5
10

ce qui me semble correct non?

[millie]

Oui, mais tu as un pas h de 2.5, c'est beaucoup trop. Il faut prendre des pas petits (au moins en dessous de 0.1). Plus h sera grand et plus ta solution divergera vite.


Par exemple, si tu fais la résolution bête :
y(0) = 1
y'(t) = y => y(t+h) - y(t) = y(t) * h => y(t+h) = y(t) (h+1)

Si h très petit (par exemple 0.01) => y(0.01) = t(0) (0.01 + 1)
Soit exp(0.01) = 1.01 (ce qui est pas loin d'être juste)

Si tu faisais directement, h = 2.5
Tu aurais exp(2.5) = 2.51 ce qui est totalement faux.

Cette méthode diverge beaucoup plus vite que celle de Runge-Kutta, mais c'est le "même principe". Au sens que si tu as un pas grand, ça part en vrille.


On peut lire sur Wikipedia (ici) que :

La méthode RK4 est une méthode d'ordre 4, ce qui signifie que l'erreur commise à chaque étape est de l'ordre de h^5, alors que l'erreur totale accumulée est de l'ordre de h^4.

L'erreur à chaque étape est de h^5. Pour h = 0.1, tu as une erreur bornée par 0.1^5, ce qui est bien. Pour ton h = 2.5, ton erreur est bornée par plus de 40...


C'est pour ça que si tu fais un module pour gérer les équations de Runge-Kutta, c'est mieux de laisser à l'utilisateur le pas de discrétisation (temporelle ou spatiale) qu'il souhaite. Plutôt que de passer par le nombre de point qu'il souhaite (car si il en donne un trop petit, il faudra de toute manière en calculer plus pour avoir quelque chose de correct).

[Razgriz]

Ouais bien vu, le truc c'est que je suisbeaucoup trop imprécis. J'ai essayé sur 1000 noeuds ça marche super bien.

Merci donc et désolé pour le dérangeemnt.

[millie]

Ouais bien vu, le truc c'est que je suisbeaucoup trop imprécis. J'ai essayé sur 1000 noeuds ça marche super bien.

Merci donc et désolé pour le dérangeemnt.

De rien, profites-en pour cliquer sur

[Razgriz]

Oups déso, j'avais oublié.
Merci pour le rappel