Error in stiff ode plot

Hi I've been having a promblem with ploting a stiff ode, I don't understand the error message that comes up as the same code has worked for a diffrent model I have investigate.

function dxdt=f(t,x) 
dxdt=zeros(5,1);
r=1; k=1; dv=1; du=1; hu=1; he=1; hv=1; delta=1; pm=1;M=1; pe=1; de=1; dt=1; omega=1; b=1; 
dxdt(1)=r*x(1)*(1-(x(1)+x(2))/k)-x(5)*dv*(1-exp(-hu*x(1)))-x(1)*du*(1-exp(-he*x(4)));
dxdt(2)=x(5)*dv*(1-exp(-hu*x(1)))-x(2)*du*(1-exp(-he*x(4)))-delta*x(2);
dxdt(3)=pm*(1-exp(-hv*x(5)))*x(3)*(1-x(3)/M);
dxdt(4)=x(3)*pe*(1-exp(hv(x(5)+x(1))))-de*x(4)-dt*x(1)*x(4);
dxdt(5)=delta*b*x(2)-omega*x(5);
end

function [T,X] = ffig()
tspan=[0 150];
x1_0=10^3; 
x2_0=0;
x4_0=0;
x5_0=1;
[T_1,X1] = ode15s(@f,tspan,[x1_0 x2_0 1 x4_0 x5_0]);

plot(T_1,X1(:,1),'k')


end

thanks in advance :)

Answer 1

[Star Strider]

You omitted an operator (that you likely intend to be a multiplication operator) ...

dxdt(4)=x(3)*pe*(1-exp(hv(x(5)+x(1))))-de*x(4)-dt*x(1)*x(4);
                         ↑ ← HERE

You also need only one ode15s call.

Try this:

function dxdt=ivl(t,x) %Xu=x(1), Xi=x(2), Xm=x(3), Xe=x(4), Xv=x(5)
dxdt=zeros(5,1);
r=0.927; k=1.8182E5; dv=0.0038E-3; du=2.0; hu=0.5E-3; he=5E-7; hv=4E-8; delta=1; pm=2.5;
M=10; pe=0.4; de=0.1; dt=5E-6; omega=2.042; b=1E6; 
dxdt(1)=r*x(1)*(1-(x(1)+x(2))/k)-x(5)*dv*(1-exp(-hu*x(1)))-x(1)*du*(1-exp(-he*x(4)));
dxdt(2)=x(5)*dv*(1-exp(-hu*x(1)))-x(2)*du*(1-exp(-he*x(4)))-delta*x(2);
dxdt(3)=pm*(1-exp(-hv*x(5)))*x(3)*(1-x(3)/M);
dxdt(4)=x(3)*pe*(1-exp(hv*(x(5)+x(1))))-de*x(4)-dt*x(1)*x(4);
dxdt(5)=delta*b*x(2)-omega*x(5);
end

tspan=[0 150];
x1_0=10^3; %per thousand 
x2_0=0;
x4_0=0;
x5_0=1;
[T,X] = ode15s(@ivl,tspan,[x1_0 x2_0 1 x4_0 x5_0]);

figure
plot(T, X)
grid

That should do what you want it to do.

[Becca Andrews]

that's grand, thanks so much for the help!

Can I just ask as now there is only one plot command, does this mean I'd need to do it several times to plot the diffrent vaules of x(3)?

[Star Strider]

As always, my pleasure!

Use a for loop:

tspan = linspace(0, 150, 50);
x3v = [1 200 300 344 345 400];
x1_0=10^3; %per thousand 
x2_0=0;
x4_0=0;
x5_0=1;
for k1 = 1:numel(x3v)
    [T,X{k1}] = ode15s(@ivl,tspan,[x1_0 x2_0 x3v(k1) x4_0 x5_0]);
end

figure
for k1 = 1:numel(x3v)
    subplot(numel(x3v), 1, k1)
    semilogy(T, X{k1})
    grid
end

Defining ‘tspan’ as I did here means that you can directly compare any or all of the solutions from any of the ‘X’ outputs with any of the others, since they all have the same times associated with them.