Discussion :

[akramakram]

Bonjour,



et merci.

[pierre.marechal]

Hello akramakram,
Voici joint un exemple de script Scilab pour résoudre une Equation Différentielle Ordinaire avec Scilab.
Il y a encore du boulot pour résoudre le problème que tu poses...

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
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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
// Nettoyage
xdel(winsid());
clear;
 
// Constantes
lw   = 3;
fs   = 4;
 
// Circuit RLC
R	= 0.2*10^(3/2);
L	= 1e-3;
C	= 1e-6;
 
// Paramètres canoniques
E	= 5;                disp(sprintf('Paramètres canoniques:'));
w0	= 1/sqrt(L*C);      disp(sprintf('   w0  = %.0f rad/s',w0));
m	= (R/2)*sqrt(C/L);  disp(sprintf('   m   = %.2f',m));
//m	= 1;                disp(sprintf('   m   = %.2f',m));
// ODE d'ordre 2 décomposée en 2 composantes d'ordre 1:
function fprime=F(t, f)
	fprime(1) = f(2);
	fprime(2) = -2*m*w0*f(2) -w0^2*f(1) + w0^2*E;
endfunction
 
// Vecteur temps
tm  = 2e-3;
nt  = 500;
t   = (0:(nt-1))/(nt-1)*tm;
// Conditions initiales
f0  = 0; //f(0)
fp0 = 0; //f'(0)
y0  = [f0;fp0]; // Vecteur colonne des conditions initiales
 
// Résolution de l'équation différentielle du 1er ordre
y  = ode(y0,0,t,F) // y est une matrice à 2 lignes
f_num  = y(1,:); // Fonction solution
fp_num = y(2,:); // Dérivée solution
 
// Solution exacte
if (0<m) & (m<1)
    f_exacte = f0 + (E-f0)*( 1 - (exp(-m*w0*t)./sqrt(1-m^2))...
                                .*cos(w0*t*sqrt(1-m^2)-atan(m/sqrt(1-m^2))));
elseif m==0
    f_exacte = f0 + (E-f0) * ( 1 - cos(w0*t) );
elseif m==1
    f_exacte = f0 + (E-f0)*( 1 - exp(-w0*t).*(1+w0*t) );
elseif m>1
    w1 = w0*(m - sqrt(m^2 - 1))
    w2 = w0*(m + sqrt(m^2 - 1))
    f_exacte = f0 + (E-f0)*( 1 - ( w2*exp(-w1*t) - w1*exp(-w2*t) )/(w2-w1) );
end
 
// Figure
cf1 = figure(1);
cf1.figure_name = 'OML2_TD3_Ex2_EDO2';
// Fond blanc
cf1.background  = color('white');
// Tracé + Grille
plot(t*1e6,E*ones(1,nt),'--k', 'LineWidth',lw+1);
e = gce(); e.children(1).line_style = 9;
plot(t*1e6,f_num,   'b', 'LineWidth',lw); xgrid();
plot(t*1e6,f_exacte,'r', 'LineWidth',lw);
e = gce(); e.children(1).line_style = 7;
 
// Axes
ca1 = get('current_axes');
ca1.font_style  = 1;
ca1.font_size   = fs;
ca1.data_bounds = [0,0;tm*1e6,E*2];
// Labels
xlabel('$ t (\mu s) $','font_size',fs+1);
ylabel('$ V_s(t) $', 'font_size',fs+1);
title('$ EDO: V_s'''' + 2.m. \omega_0.V_s'' + \omega_0^2 V_s = \omega_0^2 E; CI: V_{s,0}=0, V''_{s,0}=0 $','font_size',fs);
 
// Légende
hl=legend('$ E $', '$ f_{num}(t) $', '$ f_{exacte}(t) $');
hl.background = 8;
hl.font_size  = fs;
 
// Enregistrement
xs2svg(cf1,'EDO2_Scilab_Fig1.svg');
 
function ys0=y_sol_y0(t,E,f0,m,w0)
    if (0<m) & (m<1)
        ys0 = f0 + (E-f0)*( 1 - (exp(-m*w0*t)./sqrt(1-m^2))...
                                    .*cos(w0*t*sqrt(1-m^2)-atan(m/sqrt(1-m^2))));
    elseif m==0
        ys0 = f0 + (E-f0) * ( 1 - cos(w0*t) );
    elseif m==1
        ys0 = f0 + (E-f0)*( 1 - exp(-w0*t).*(1+w0*t) );
    elseif m>1
        w1 = w0*(m - sqrt(m^2 - 1));
        w2 = w0*(m + sqrt(m^2 - 1));
        ys0 = f0 + (E-f0)*( 1 - ( w2*exp(-w1*t) - w1*exp(-w2*t) )/(w2-w1) );
    end
endfunction         
 
// Conditions Initiales testées
//f0k = [-10 -6 -2 +2 +6 +10];
f0k = (-2:+2)*E;
nf0 = length(f0k);
 
// Figure
cf2 = figure(2);
cf2.figure_name = 'OML2_TD3_Ex1_EDO2_f0';
cf2.background  = 8;
// Colormap
nc	= 255;
cm	= jetcolormap(nc);
 
// Tracé et légende
legstr = '';
for k=1:nf0
    fk(k,:) = y_sol_y0(t,E,f0k(k),m,w0);
    ind = fix(nc*k/nf0);
    plot(t*1e6,fk(k,:),'Color',cm(ind,:),'LineWidth',lw);
    legstr=sprintf('%s;""$ V_{s,0} = %d $""',legstr,f0k(k));
end
xgrid;
// Axes
ca2 = get('current_axes');
ca2.font_style  = 1; //8ca1.font_style  = 1;
ca2.font_size   = fs;
 
// Labels
xlabel('$ t (\mu s) $','font_size',fs+1);
ylabel('$ V_s(t) $', 'font_size',fs+1);
title('$ EDO: V_s'''' + 2. \omega_0.V_s'' + \omega_0^2 V_s = \omega_0^2 E; CI: V''_{s,0}=0 $','font_size',fs);
// Légende
legs = part(legstr, 2:length(legstr));
hl=legend(evstr(legs));
hl.background = 8;
hl.font_size  = fs;
hl.legend_location = 'in_lower_right';
 
// Enregistrement
xs2svg(cf2,'EDO2_Scilab_Fig2_f0.svg');

Tu nous feras un retour si tu y arrives ?