Error in stiff ode plot

10 次查看（过去 30 天）
Becca Andrews
Becca Andrews 2019-3-16
编辑： Becca Andrews 2019-3-24
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 :)

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2019-3-16
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 个评论
Becca Andrews
Becca Andrews 2019-3-16
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
Star Strider 2019-3-16
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by