Discussion :

[osimte]

Bonjour,
Je suis en train de résoudre (ou bien, essayer de..) ce système


a(t) est une fonction gaussienne qui décrit les ‘inputs’ du système


J’ai fait comme ça:

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
function dx=f(t, x, a)
    dx (1) = a(t) - Pgamma*x(2)*(1-R) // dg/dt
    dx (2) = x(3) - Pbeta*x(2) // di/dt
    dx (3) = Pdelta*x(1) - Peta*x(1)^2 // dq/dt
endfunction
function [res]=a(t)
    inter = 2; 
    larg = 1; 
    haut = 220; 
    res = zeros(1,length(t))
    for i=1:length(t)
 
        if modulo(t(i),inter)<=larg/2 | modulo(t(i),inter)>=inter-larg/2 then
            if modulo(t(i),inter)>=inter-larg/2 then
                    x = modulo(t(i),inter) - (inter-larg/2);
            elseif modulo(t(i),inter)<=larg/2
                    x = modulo(t(i),inter) + larg/2;
            end
            res(i) = haut*exp(-(x-larg/2)^2/(larg/(24)));
        else
            res(i) = 0;
        end
    end
endfunction
Pgamma = 1.82; 
Pbeta = 1/25; 
Peta = 1.05; 
Pdelta = 90*Peta; 
R = 0;
X0=[90;25;1]; //[g0;i0;q0]
t0=0;
t=0:0.005:3; 
x=ode(X0,t0,t,f)

Mais je reçois une erreur, qui disparait si j’enlève le terme au carré de la troisième équation:

lsoda--  at t (=r1), mxstep (=i1) steps   
needed before reaching tout
      where i1 is :        500                                                  
      where r1 is :   0.1132577023323D+01      

Est-ce-que quelqu’un a des suggestions ?

Merci beaucoup

[phryte]

Bonjour,
Une idée :

Code :

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

 dx (3) = Pdelta*x(1) - Peta*x(1)*x(1)) // dq/dt

[osimte]

Bonjour,
Une idée :

Code :

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

 dx (3) = Pdelta*x(1) - Peta*x(1)*x(1)) // dq/dt

Bonjour,
J'avais essayé mais ça ne change rien. J'ai aussi essayé d'utiliser une fonction externe, mais ça ne change rien non plus...

Par contre j'ai noté que ce n'est pas un problème lié au carré: en enlevant, par exemple, le dernier terme dans la première équation, j'ai un résultat (pas correct, mais je n'ai pas d'erreur).

[phryte]

Limite de fonctionnement :

Code :

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

t=0:0.005:1.1;

[osimte]

Limite de fonctionnement :

Code :

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

t=0:0.005:1.1;

Pourquoi? Ça vient de où?
Excuse-moi pour la question naïf...

[phryte]

Bonjour,
Le problème vient de l'initialisation du maximum d'itérations des calculs à 500 :

Code :

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

odeoptions = [1,0,0,%inf,0,2,500,12,5,0,-1,-1]

Je n'ai pas réussi à forcer davantage.

[rouxph]

Le problème ne vient pas de scilab mais du problème considéré. Dans la solution la 3ième fonction "explose" au-delà de t=1 , trace le graphe pour vérifier :

Code :

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

1
2
3
4
// je ne modifie que les dernières lignes du code
t=0:0.005:1; // limiter t<=1
x=ode(X0,t0,t,f);// résoudre
plot(t,x(1,:),'-r',t,x(2,:),'-b',t,x(3,:),'-g')// afficher les graphes

la courbe en vert atteint déjà des valeurs de l'ordre de 10^7 pour t=1