Résolution d'équation différentielle RK4

Version imprimable

Bonjour à tous,

Je sollicite votre aide car je suis complètement perdue et une vraie débutante!

J'ai un projet en C++ et voici mon sujet:

On considère un câble sous tension auquel sont rigidement et régulièrement attachés des pendules. Les pendules sont couplés grâce au câble (à travers sa constante de torsion). Dans un tel système on peut observer une large gamme de phénomènes ondulatoires. Le but de cet projet est d’étudier une solution très particulière : le soliton.

Imaginons qu’une des extrémités du câble est attachée à une manivelle qui peut tourner librement. Il est alors possible de donner une impulsion au système en faisant un tour rapide ce qui déclenche la propagation d’un soliton. Dans ce projet, on considérera les pendules individuellement. Il n’est pas demandé de passer au modèle continu et de résoudre l’équation obtenue. Pour chaque pendule n, l’équation d’évolution s’écrit : d²(theta)/dt²=w0²sin(theta)-beta[theta(i+1)+theta(i-1)-theta(i)]

On résoudra cette équation pour chaque pendule. En donnant un « tour de manivelle numérique », on essayera d’obtenir la solution soliton
Code:

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 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 #include <iostream> #include <fstream> #include <cmath> using namespace std; double f(double y) { return y; } double g( double theta[], int N,int i, double pas, double w0, double beta) { return -w0*w0*sin(theta[i]+pas)-beta*(theta[i+1]+theta[i-1]-theta[i]); } void RK4(double t,double theta[], double vit[], int N) { double h=0.01, tmp; double w0=1.0; double beta= 1.0; for (int i=0; i<N; i++) { double pas=0.0; double k1=g(theta,N,i,pas,w0,beta); pas= (h/2.)*k1; double k2=g(theta, N,i,pas, w0,beta); pas= (h/2.)*k2; double k3=g(theta,N,i,pas, w0,beta); pas= h*k3; double k4=g(theta, N,i,pas, w0,beta); double a=(k1+2*k2+2*k3+k4)/6.; double l1=f(vit[i]); double l2=f(vit[i]+(h/2.)*l1); double l3=f(vit[i]+(h/2.)*l2); double l4=f(vit[i]+ h*l3); double b=(l1+2*l2+2*l3+l4)/6.; theta[i]= theta[i]+h*a; vit[i] = vit[i] + h*b; } } int main () { int N,i; double t=0.,tmp=0,h=0.01; cout << "Nombre de pendules N = "; cin >> N; double theta[N], vit[N]; theta[0]=0.01; vit [0]=0; //cout << "entrer la valeur de h : "; // cin >> h; fstream w; w.open("soliton.txt",ios::out); double w0=1.0; double beta= 1.0; do { RK4(t, theta, vit ,N); w<<t << " "<< theta[i] << " "<< vit[i] << endl; } t=t+h; } while (t<2.); w.close(); return 0; }
J'ai écris ce programme jusqu'à présent (je sais bien qu'il y a énormément de soucis avec mais je ne sais pas comment faire).

J'espère que vous allez m'aider au moins à me diriger vers quelque chose qui marche!

Merci beaucoup!!:roll:

Bonjour,
Je n'ai rien compris au problème physique. Ce qui est sûr c'est que ton code est pour le moins incomplet.
On a une équation qui décrit le mouvement de chaque pendule en fonction de paramètres.
Il faut tout d'abord bien voir à quoi correspond chacun des paramètres de l'équation.
* i c'est quoi? Si c'est l'indice d'un pendule, on aura peut-être besoin du nombre N de pendules.
* w0 c'est quoi?
* theta est surement un angle et theta(i) est-il une fonction de i ou un indexage theta_i?
* t doit être le temps, et certainement que le programme doit fournir un résultat pour divers quanta de temps
* n'étant pas mathématicien, je ne sais pas quoi faire de la dérivée seconde.
Ensuite il faut certainement propager et itérer sur les valeurs que l'on connait.

Pour une fonction plus simple, par exemple dx/dt = i / (x+1) je coderais :
Code:

1 2 3 double fctX( double dt , int i , double prevX ) { return prevX + ( i / (prevX+1) ) * dt; }
Et ensuite je ferai une double boucle (sur dt et sur i) pour parcourir les cas possibles et calculer les valeurs à chaque instant :
Code:

1 2 3 4 5 6 7 8 9 10 11 int main() { std::array<double,10> lesX{}; // on met les 10 positions initiales à 0 double dt = 0.1; // incrémentation du temps for ( double t = dt ; t < 1. ; t += dt ) { std::cout << "À l'instant de date " << t << "secondes\n"; for ( int i = 0 ; i < lesX.size() ; ++i ) { // pour chacun des éléments lesX[i] = fctX( dt , i , lesX[i] ); // on calcule sa position std::cout << "la position x du " << (i+1) << "ème élément est " << lesX[i] << " mètres\n"; } } }

28/03/2016, 16h54
Whira

Bonjour,
Merci pour ta réponse.

Alors i est bien l'indice de la position du pendule. i est une position et i-1 la position d'avant et i+1 la position d'après. Effectivement, il y a bien N, nombre total des pendules, qui doit entrer en jeu et c'est là que je bloque aussi.
w0 et beta sont juste des paramètres avec w0=-mgl/I et avec m masse des pendules, g gravité, l longueur des pendules et I moment d'inertie du pendule.
theta est effectivement l'angle du pendule et le i est un indexage de position.

La fonction à appeler devrait avoir la forme suivante (attention j'ai posé dTheta/dt = ... au lieu de d²theta/d²t, ce qui me permet d'incrémenter directement theta_i)
Code:

1 2 3 4 5 6 7 8 // la table theta est initialisée avec l'angle initial, theta[1] à theta[n] évoluent et theta[n+1] mis à 0 // i est valide de 1 à n pour la fonction double fct( double theta[], int i, double dt, double w0, double beta ) { double delta = -w0*w0*sin(theta[i]) - beta * (theta[i+1] +theta[i-1] -theta[i]); theta[i] += delta * dt; // surement faux! return theta[i]; }
Après il faut voir quelles sont les itérations d'appel à effectuer, et l'effet d'être en d²theta/d²t

29/03/2016, 10h57
Whira

La boucle sur le i pour les N pendules se fera dans le main alors?

Si je comprend bien, tu veux calculer pour chaque pendule les valeurs kx pour ensuite en faire une approximation de theta.

Dans un premier temps, je te déconseille d'utiliser les tableaux à la C dans un programme C++. Utilise plutôt les conteneurs offerts par la STL à savoir std::array (pour des tableaux de taille fixe) ou std::vector.

Si je reprend ton équation différentielle :

Code:

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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include <iostream>
#include <array>
#include <vector>
#include <cmath>
 
double g(const std::vector<double> &theta, int i, double pas, double w0, double beta)
{
    return -w0*w0*sin(theta[i] + pas) - beta*(theta[i+1] + theta[i-1] + theta[i]);
}
 
int main()
{
    // constantes
    constexpr double W0 = 1.0;
    constexpr double BETA = 1.0;
    constexpr double DT = 0.01;
    constexpr double H = 0.01;
 
    // pentes de RK4
    std::array<double, 4> k;
 
    // nombre de pendules
    size_t n = 0;
    do
    {
        std::cout << "Nombre de pendules N = ";
        std::cin >> n;
    }
    while (n < 1);
    std::vector<double> theta(n, 0.);
 
    // angle initial
    theta[0] = 0.01;
 
    // pour chaque variation du temps
    for (double t = DT; t < 2.; t += DT)
    {
        // pour chaque pendule
        for (size_t i = 0; i < n; i++)
        {
            // calcul des k
            k[0] = g(theta, i, 0., W0, BETA);
            k[1] = g(theta, i, (H/2.)*k[0], W0, BETA);
            k[2] = g(theta, i, (H/2.)*k[1], W0, BETA);
            k[3] = g(theta, i, H*k[2], W0, BETA);
 
            // calcul de l'angle approximé
            theta[i] += H*(k[0] + 2.*k[1] + 2.*k[2] + k[3])/6.0;
            std::cout << "theta[" << i << "](" << t << ") = " << theta[i] << std::endl;
        }
    }
 
    return 0;
}

Que se passe-t-il pour i = 0 et i = n-1 ?

29/03/2016, 11h26
Whira

Je te remercie pour ton aide.

Le programme que je voudrais pouvoir faire doit être capable de calculer lui même les angles theta de chaque pendule i à chaque instant t à partir d'une position initiale de l'angle. Le but serait aussi de tracer les solutions et de les comparer avec la solution théorique attendue.
29/03/2016, 11h39
RPGamer

C'est ce que j'avais compris mais tu ne réponds pas à ma question ;)
29/03/2016, 12h00
Whira

Position initiale au repos sinon ... je ne sais pas trop justement!

Dans mon exemple, remplace

Code:

1
2
3
4
double g(const std::vector<double> &theta, int i, double pas, double w0, double beta)
{
    return -w0*w0*sin(theta[i] + pas) - beta*(theta[i+1] + theta[i-1] + theta[i]);
}

par

Code:

1
2
3
4
double g(const std::vector<double> &theta, int i, double pas, double w0, double beta)
{
    return -w0*w0*sin(theta.at(i) + pas) - beta*(theta.at(i+1) + theta.at(i-1) + theta.at(i));
}

et observe ce qui se passe à l'exécution.