Six coupled differential equations

New to matlab & trying to solve these equations for y(1), y(2), etc. values. I have six unknowns and six equations - why would this script fail?
function dY = rate( t,Y )
%
dY = zeros(6,1);
A51 = 26993 ;
A52 = 9894 ;
A53 = 897 ;
A54 = 1729 ;
A41 = 455 ;
A42 = 147;
A31 = 398;
A65 = 392937;
A43 = 102819;
T3 = 240E-6;
A32 = 1/T3 - A31;
C22= 10^-198;
C33= 10^-18;
R31 = 980E-9;
R13 = 1E-20;
A21 = 3094;
R36 = 240E-6;
n0=1E-20;
dY(1) = -R13*Y(1) +R31*Y(3) +A21*Y(2) +A31*Y(3) + A41*Y(4) +A51*Y(5) +C22*Y((2))^2 +C33*Y((3))^2;
dY(2) = -A21*Y(2) +A32*Y(3) +A42*Y(4) +A52*Y(5) -2*C22*(Y(2))^2;
dY(3) = R13*Y(6) -R31*Y(3) -R36*Y(3) - A31*Y(41) - A32*Y(3) + A43*Y(4) + A53*Y(5) -9*C33*(Y(3))^2;
dY(4) = -A41*Y(2) - A42*Y(4) - A43*Y(4) + A54*Y(2) +C22*(Y(2))^2;
dY(5) = -A51*Y(5) - A52*Y(5) - A53*Y(5) - A54*Y(5) + A65*Y(6);
dY(6) = R36*Y(3) - A65*Y(6) + C33*(Y(3))^2;
n0=Y(1)+Y(2)+Y(3)+Y(4)+Y(5)+Y(6);
end
* _MAIN_*
TSPAN = [0 0 0 0 0 0 0];
t0 = 0; tf = 20E-2;
[T,Y] = ode15('rate_eq2',TSPAN,t0,tf);
subplot(6,1,1)
plot(t,y(:,6),t,y(:,5),t,y(:,4),t,y(:,3),t,y(:,2),t,y(:,1))
legend('6','5','4','3','2','1')
title('Population different levels')

 Respuesta aceptada

Star Strider
Star Strider el 26 de Oct. de 2016
Your function is not the problem. You’re not calling it correctly, and you’re not using the correct solver.
Call it this way instead, and it works:
icv = ones(6,1); % <== CHANGE TO APPROPRIATE INITIAL CONDITIONS
t0 = 0; tf = 2E-2;
TSPAN = linspace(t0, tf, 150);
[T,Y] = ode15s(@rate_eq2,TSPAN,icv);
figure(1)
plot(T, Y)
grid
legend('6','5','4','3','2','1', 'Location','E')
title('Population different levels')
The ‘icv’ assignment are the initial conditions vector. If they’re all initially zero, then the solution is zero, so at least some of them have to be non-zero. You were also confusing the initial conditions and the time span arguments.

2 comentarios

Grace
Grace el 26 de Oct. de 2016
Great, thank you so much for the help!
Star Strider
Star Strider el 26 de Oct. de 2016
My pleasure!

Iniciar sesión para comentar.

Más respuestas (0)

Etiquetas

Preguntada:

el 26 de Oct. de 2016

Comentada:

el 26 de Oct. de 2016

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by