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 :)
采纳的回答
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
2019-3-16
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 个）
类别
在 帮助中心 和 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 
