Discussion :

[elmoukrie]

Bonjour,

Je souhaite utiliser un logiciel professionnel sur lequel on peut charger un code MATLAB. Ce dernier doit être rédigé d'une façon "spéciale".

Sur un tutorial descriptif, on met un exemple d'un système de deux ODE (qui est joint) avec le code MATLAB correspondant qu'il faut charger dans le logiciel.

Mon cas concerne 11 ODE, j'ai refait le code en suivant l'exemple donné à la lettre et il est en bas. Mais il me faut quelques dernières retouches pour que le code soit valide. J'ai mis mes questions en commentaire %.

Il sera facile pour quelqu'un ayant la logique de programmation (surtout en vectoriel) sur MATLAB de voir ce qui manque ou ce qui ne va pas.


Pourriez vous SVP jeter un coup d’œil sur l’exemple et vous en inspirer pour compléter le code de ma fonction?

Mille merci à vous.

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
 
 
function [y_out varargout] = user_defined(x_in, varargin)
n = size(x_in, 1); 
N = size(x_in, 2); 
t_step = 0:10:10000; 
m1 = 11; 
m2 = length(t_step);
varargout{1} = m1; 
varargout{2} = m2;
 
p = x_in;
 
y0(1,:) = 0*ones(1, N);    % IS IT CORRECT TO PUT IT IN THIS WAY = 0*   ?
y0(2,:) = CHP_initial*ones(1, N);    % The start value of y2 is the sixth input parameter. HOW To PUT IT?
y0(3,:) = 0*ones(1, N); 
y0(4,:) = ((1-(0.14*CHP_initial)*0.24285)*1000)/18*ones(1, N);    % The start value of y4 is A FUNCTION OF the sixth input parameter: ((1-(0.14*CHP_initial*0.24285)*1000)/18.  Is it correct like this?
y0(5,:) = 0.5*ones(1, N); 
y0(7,:) = 0*ones(1, N); 
y0(8,:) = 0*ones(1, N);
y0(9,:) = 0*ones(1, N); 
y0(10,:) = Tj0*ones(1, N);   % The start value of y10 is the fifth input parameter. HOW To PUT IT?
y0(11,:) = 0.26*ones(1, N); 
 
y0 = reshape(y0, m1*N, 1);
[t y] = ode23s(@semibatch, t_step, y0, [], p);
y_1 = reshape(y, m2, m1, N);     % What does it mean? should we keep it?
y_out = permute(y_1, [2,1,3]);   % What does it mean? should we keep it?
 
 
 
function dydt = semibatch(tin, yin, pin)
dydt = zeros(2*size(pin,2), 1);
 
 
if t > tadd     % HOW AND WHERE should I put the If test?
       F=0;
end
 
j = 1;
 
 
for i = 1:size(pin,2)
 
    dydt(j)=[(1/(1+((1-((y(j+10)-0.12)/y(j+10)))/(((y(j+10)-0.12)/y(j+10))*9))))*((pin(1,i)/(y(j+10)-0.12))*(24-y(j))-((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))+((0.0009*exp(-(2429.636)*((1/y(j+9))-(0.0029411))))*y(j+2))+(1-((y(j+10)-0.12)/y(j+10)))*(((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6))-((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3)))))/((y(j+10)-0.12)/y(j+10)));
        dydt(j+1)= -((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))-(pin(1,i)*y(j+1)/(y(j+10)-0.12));
        dydt(j+2)= ((-pin(1,i)*y(j+2)/(y(j+10)-0.12))+((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))-((0.0009*exp(-(2429.636)*((1/y(j+9))-(0.0029411))))*y(j+2))-((0.001*exp(-(2405.581)*((1/y(j+9))-(0.0029411))))*y(j+2))-(1-((y(j+10)-0.12)/y(j+10)))*(((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6))+((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3)))))/((y(j+10)-0.12)/y(j+10)));
        dydt(j+3)= ((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))+((0.001*exp(-(2405.581)*((1/y(j+9))-(0.0029411))))*y(j+2))-(PIN(1,I)*y(j+3)/(y(j+10)-0.12))-(1-((y(j+10)-0.12)/y(j+10)))*((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))/((y(j+10)-0.12)/y(j+10)); 
        dydt(j+4)= -((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4));
        dydt(j+5)= -((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5)); 
        dydt(j+6)= ((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))-((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6)); 
        dydt(j+7)= (((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6)))-((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))-((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3))))-((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3))));
        dydt(j+8)= ((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))+((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3))))+((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*(0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3)))); 
        dydt(j+9)=(1/(((y(j+10)-0.12)*1.00+0.12*0.93)*2000))*((-(y(j+10)-0.12)*(((0.15*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*(sqrt(y(j)/y(j+3)))*(y(j)*y(j+1)-(1/(0.96*exp((671.157)*((1/y(j+9))-(0.003298)))))*y(j+2)*y(j+3)))*-5580+((0.001*exp(-(2405.581)*((1/y(j+9))-(0.0029411))))*y(j+2))*-359000+((0.0009*exp(-(2429.636)*((1/y(j+9))-(0.0029411))))*y(j+2))*-163000)-0.12*((((0.00576*exp(-(7409.189)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+4))+((0.00437*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+5))+((0.004*exp(-(1804.1857)*((1/y(j+9))-(0.0029411))))*y(j+2)*y(j+6)))*-230000+(((0.000237*exp(-(18041.857)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*(sqrt(y(j)*y(j+3))))+((0.00339*exp(-(5063.74789)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j)*(sqrt(y(j)/y(j+3))))+((0.0592*exp(-(8419.53331)*((1/y(j+9))-(0.0029411))))*( 0.0017029)*y(j+7)*y(j+2)*(sqrt(y(j)/y(j+3)))))*-90000))+pin(3,i)*(pin(5,i)-y(j+9))+24*pin(1,i)*20*(pin(4,i)-y(j+9))); 
        dydt(j+10)= pin(1,i);
        j=j+11;
 
        end