Matlab - Error en índices sistema EDOs

   
Vista:
Imágen de perfil de Geroge

Error en índices sistema EDOs

Publicado por Geroge (1 intervención) el 14/07/2014 01:44:40
Hola, pretendo resolver un sistema de EDO´s pero no he podido. Les dejo mi código:

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
function dy = concentraciones(t,y)
 
% Constantes
VI=2*10^(-5);
VII=11*10^(-5);
VIII=8*10^(-4);
Def=1.7*10^(-10);
S1=2.83*10^(-11);
N2=11944.5;
S2=2.37*10^(-5);
RfI=1.5*10^(-6);
RfII=1.375*10^(-3);
RfIII=RfI+RfII;
Keq1=106;
k1=2.5*10^(-5);
Keq2=10;
k2=1.1*10^(-6);
kA=3.0*10^(-5);
V=(VI+VII)/VII;
 
% Dominio de resolución
dr=0.5*RfI;
dt=10^(-3);
RI=dr:dr:RfI;
RII=RfI+dr:dr:RfII;
RIII=RfII+dr:dr:RfIII;
R=[RI RII RIII];
r1=length(RI);
r2=length(RI)+length(RII);
r3=length(R);
 
 
 
dy=ones(length(R),1);%matriz de concentraciones
 
 
for r=1:1:r3
 if r<r1
     v2=k2*(y(r1+1)*y(r1)-(Keq2*y(r1-1)*(0.055-2*y(r1+1))^2)/y(r1));
     dy(r)=(S1/VI)*v2;
 
 else
     if r==r1
         y(r)=2*7.65*10^(-4)+2.5-2*y(r-1);
     else
         if r==r1+1
             v2=k2*(y(r)*y(r1)-(Keq2*y(r1-1)*(0.055-2*y(r))^2)/y(r1));
             dy(r)=(V*Def/dr^2)*(2*y(r+1)-2*y(r))-(S1/VII)*v2;
 
         else
             if (r1+1<r)&&(r<r2)
                 v2=k2*(y(r1+1)*y(r1)-(Keq2*y(r1-1)*(0.055-2*y(r1+1))^2)/y(r1));
                 dy(r)=(V*Def/dr^2)*(y(r+1)-2*y(r)+y(r-1))-(S1/VII)*v2;
 
             else
                 if r==r2
                     v2=k2*(y(r1+1)*y(r1)-(Keq2*y(r1-1)*(0.055-2*y(r1+1))^2)/y(r1));
                     dy(r)=(2*V*Def/dr^2)*((dr*kA/Def)*(y(r2+2)-y(r2+1))-y(r)+y(r-1))-(S1/VII)*v2;
 
                 else
                     if r==r2+1
                         v1=k1*((y(r)*(0.055-2*y(r2))/(10^(-3)))-((y(r2)*10^(-3))/(Keq1*(0.055-2*y(r2)))));
                         y(r)=y(r+1)+((N2*S2)/(VIII*kA))*v1;
 
                     else
                         if r2+1<r
                             dy(r)=-(N2*S2/VIII)*kA*(y(r)-y(r2+1));
                         end
                     end
                 end
             end
         end
     end
 end
end
 
 
end
 
%Para resolver utilizo:
 
RfI=1.5*10^(-6);
RfII=1.375*10^(-3);
RfIII=RfI+RfII;
dr=0.5*RfI;
RI=dr:dr:RfI;
RII=RfI+dr:dr:RfII;
RIII=RfII+dr:dr:RfIII;
R=[RI RII RIII];
 
y0=zeros(1,length(R));
for i=1:length(R)
    if i==1
        y0=0;
    else
        if i==2
            y0=2.5;
        else
            if i==length(R)-1
                y0=0;
            else
                if i==length(R)
                    y0=7.65*10^(-4);
                end
            end
        end
    end
end
 
 
[t,y] = ode45(@concentraciones,[0 50],y0);

El error que arroja matlab es:

Attempted to access y(3); index out of bounds because numel(y)=1.

Error in concentraciones (line 40)
v2=k2*(y(r1+1)*y(r1)-(Keq2*y(r1-1)*(0.055-2*y(r1+1))^2)/y(r1));

Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in Solve_concentraciones (line 30)
[t,y] = ode45(@concentraciones,[0 50],y0);

No entiendo el error. Alguna sugerencia por favor?
Muchas gracias
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder