Error in stiff ode plot
1 view (last 30 days)
Show older comments
Becca Andrews
on 16 Mar 2019
Edited: Becca Andrews
on 24 Mar 2019
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 :)
0 Comments
Accepted Answer
Star Strider
on 16 Mar 2019
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.
2 Comments
Star Strider
on 16 Mar 2019
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.
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!