Borrar filtros
Borrar filtros

Error while solving a coupled system of ODE

1 visualización (últimos 30 días)
F.O
F.O el 17 de Nov. de 2019
Comentada: Star Strider el 17 de Nov. de 2019
I am trying to solve a system of 6 coupled ODE my function works but I dont know why calling ode45 doesnt work.
My function is:
function [dxdt] = DM_bacteria(t,x,sigma,r,f,s,epsilon,omega,gamma12,g12)
dxdt(4)=x(4)*(r(4)*(1-(sigma(4,4)*x(4)+sigma(4,1)*x(1)+sigma(4,2)*x(2)+sigma(4,3)*x(3))/r(4))-f(4,1)*x(5)-f(4,2)*x(6))+s(4)
dxdt(5)=x(5)*(-epsilon(1)*(1+(omega(1,1)*x(5)+omega(1,2)*x(6))/epsilon(1))-g12*x(5)+beta(1,1)*f(1,1)*x(1)+beta(2,1)*f(2,1)*x(2)+beta(3,1)*f(3,1)*x(3)+beta(4,1)*f(4,1)*x(4))
dxdt(6)=x(6)*(-epsilon(2)*(1+(omega(2,2)*x(6)+omega(2,1)*x(5))/epsilon(2))-gamma12*g12*x(5)+beta(1,2)*f(1,2)*x(1)+beta(2,2)*beta(2,1)*f(2,2)*x(2)+beta(3,2)*f(3,2)*x(3)+beta(4,2)*f(4,2)*x(4))
dxdt(1)=x(1)*(r(1)*(1-(sigma(1,1)*x(1)+sigma(1,2)*x(2)+sigma(1,3)*x(3)+sigma(1,4)*x(4))/r(1))-f(1,1)*x(5)-f(1,2)*x(6))+s(1)
dxdt(2)=x(2)*(r(2)*(1-(sigma(2,2)*x(2)+sigma(2,1)*x(1)+sigma(2,3)*x(3)+sigma(2,4)*x(4))/r(2))-f(2,1)*x(5)-f(2,2)*x(6))+s(2)
dxdt(3)=x(3)*(r(3)*(1-(sigma(3,3)*x(3)+sigma(3,1)*x(1)+sigma(3,2)*x(2)+sigma(3,4)*x(4))/r(3))-f(3,1)*x(5)-f(3,2)*x(6))+s(3)
end
The script :
sigma=ones(4,4)
f=ones(4,2)
g12=1
gamma12=2
beta=ones(4,2)
omega=ones(2,2)
epsilon=[4,2]
s=[0 0 0 0]
x0=[600 600 600 600 600 600]
t = [0 20];
r=[1 2 3 4]
t=linspace(0,20,1000)
[dxdt] = DM_bacteria(t,x0,sigma,r,f,s,epsilon,omega,gamma12,g12)
[t,x]=ode45(@DM_bacteria,t,x0)
I got error :(I have all the argument in the function)!!
Not enough input arguments.
Error in DM_bacteria (line 3)
dxdt(4)=x(4)*(r(4)*(1-(sigma(4,4)*x(4)+sigma(4,1)*x(1)+sigma(4,2)*x(2)+sigma(4,3)*x(3))/r(4))-f(4,1)*x(5)-f(4,2)*x(6))+s(4)
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in DynamicalModelBactera (line 14)
[t,x]=ode45(@DM_bacteria,t,x0)

Respuesta aceptada

Star Strider
Star Strider el 17 de Nov. de 2019
Your ‘DM_bacteria’ function must return a column vector output.
This worked when I ran it just now:
function [dxdt] = DM_bacteria(t,x,sigma,r,f,s,epsilon,omega,gamma12,g12)
dxdt(4,:)=x(4)*(r(4)*(1-(sigma(4,4)*x(4)+sigma(4,1)*x(1)+sigma(4,2)*x(2)+sigma(4,3)*x(3))/r(4))-f(4,1)*x(5)-f(4,2)*x(6))+s(4);
dxdt(5,:)=x(5)*(-epsilon(1)*(1+(omega(1,1)*x(5)+omega(1,2)*x(6))/epsilon(1))-g12*x(5)+beta(1,1)*f(1,1)*x(1)+beta(2,1)*f(2,1)*x(2)+beta(3,1)*f(3,1)*x(3)+beta(4,1)*f(4,1)*x(4));
dxdt(6,:)=x(6)*(-epsilon(2)*(1+(omega(2,2)*x(6)+omega(2,1)*x(5))/epsilon(2))-gamma12*g12*x(5)+beta(1,2)*f(1,2)*x(1)+beta(2,2)*beta(2,1)*f(2,2)*x(2)+beta(3,2)*f(3,2)*x(3)+beta(4,2)*f(4,2)*x(4)) ;
dxdt(1,:)=x(1)*(r(1)*(1-(sigma(1,1)*x(1)+sigma(1,2)*x(2)+sigma(1,3)*x(3)+sigma(1,4)*x(4))/r(1))-f(1,1)*x(5)-f(1,2)*x(6))+s(1);
dxdt(2,:)=x(2)*(r(2)*(1-(sigma(2,2)*x(2)+sigma(2,1)*x(1)+sigma(2,3)*x(3)+sigma(2,4)*x(4))/r(2))-f(2,1)*x(5)-f(2,2)*x(6))+s(2);
dxdt(3,:)=x(3)*(r(3)*(1-(sigma(3,3)*x(3)+sigma(3,1)*x(1)+sigma(3,2)*x(2)+sigma(3,4)*x(4))/r(3))-f(3,1)*x(5)-f(3,2)*x(6))+s(3);
end
sigma=ones(4,4);
f=ones(4,2);
g12=1;
gamma12=2;
beta=ones(4,2);
omega=ones(2,2);
epsilon=[4,2];
s=[0 0 0 0];
x0=[600 600 600 600 600 600];
t = [0 20];
r=[1 2 3 4];
t=linspace(0,20,1000);
% [dxdt] = DM_bacteria(t,x0,sigma,r,f,s,epsilon,omega,gamma12,g12)
[t,x]=ode45(@(t,x)DM_bacteria(t,x,sigma,r,f,s,epsilon,omega,gamma12,g12),t,x0);
figure
plot(t, x)
grid
ylim([-10, 10])
Experiment to get different results.
  2 comentarios
F.O
F.O el 17 de Nov. de 2019
Should I use ode45 or ode15s here ?
Star Strider
Star Strider el 17 de Nov. de 2019
Your system does not appear to be ‘stiff’, and ode45 had no problems with it. I would just continue with ode45, unless the parameters change significantly in their relative magnitudes and the system becomes ‘stiff’.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by