[MATLAB数学相关] ode求解微分方程

[chenzhongsheng]

求助：怎么我设置的区间是[0 500]，但是画出图的横坐标不是0到500，还有怎么把odefun里面的中间变量传递出来？？我的程序是这样的：

%在second.m里面
clear all;
clc


s0=struct('y',{'ax1';'bx1';'am';'bm';'ah';'bh';'aj';'bj';'ad';'bd';'af';'bf'}...
,'c1',{'p23';'p27'; '0'; 'p34';'p36';'p38';'p41';'p45';'p48';'p52';'p56';'p60'}...
,'c2',{'p24' ;'p28' ;'0' ;'p35';'p37';'0';'p42';'0';'p49';'p53';'p57';'p61'}...
,'c3',{'0';'0';'p31';'0';'0';'0';'0';'0';'0';'0';'0';'0';}...
,'c4',{'0';'0';'p32';'0';'0';'0';'0';'0';'0';'0';'0';'0';}...
,'c5',{'p25';'p29';'p3';'0';'0';'p39';'p43';'p46';'p50';'p54';'p58';'p62'}...
,'c6',{'p26';'p30';'-1';'0';'0';'p40';'p44';'p47';'p51';'p55';'p59';'p63'});

s=set_struct(s0);

options=odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6],...
'outputfcn','odephas2');


[T,Y]=ode23t(@deal_diff,[0 500],[1 2 3 4 5 6 7 8 ]',options,s,0.000002);


plot(T,Y(:,7),T,Y(:,8))



%set_struct.m里面
function a = set_struct(a)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here

a(1).c1=0.0018;
a(2).c1=8.7142e-4;
a(3).c1=0;
a(4).c1=0.7096;
a(5).c1=5.4980e-10;
a(6).c1=10.7578;
a(7).c1=0.0011;
a(8).c1=7.3598;
a(9).c1=0.0697;
a(10).c1=0.0037;
a(11).c1=1.4383e-4;
a(12).c1=1.4391;



a(1).c2=0.083;
a(2).c2=-0.06;
a(3).c2=0;
a(4).c2=-0.056;
a(5).c2=-0.25;
a(6).c2=0;
a(7).c2=-0.25;
a(8).c2=0;
a(9).c2=-0.01;
a(10).c2=-0.017;
a(11).c2=-0.008;
a(12).c2=-0.02;


a(1).c3=0;
a(2).c3=0;
a(3).c3=-1;
a(4).c3=0;
a(5).c3=0;
a(6).c3=0;
a(7).c3=0;
a(8).c3=0;
a(9).c3=0;
a(10).c3=0;
a(11).c3=0;
a(12).c3=0;


a(1).c4=0;
a(2).c4=0;
a(3).c4=47;
a(4).c4=0;
a(5).c4=0;
a(6).c4=0;
a(7).c4=0;
a(8).c4=0;
a(9).c4=0;
a(10).c4=0;
a(11).c4=0;
a(12).c4=0;



a(1).c5=0.057;
a(2).c5=-0.04;
a(3).c5=85;
a(4).c5=0;
a(5).c5=0;
a(6).c5=-0.082;
a(7).c5=-0.2;
a(8).c5=-0.1;
a(9).c5=-0.072;
a(10).c5=0.05;
a(11).c5=0.15;
a(12).c5=-0.2;



a(1).c6=0.0578;
a(2).c6=2.2255;
a(3).c6=-1;
a(4).c6=0;
a(5).c6=0;
a(6).c6=6.3281;
a(7).c6=5.9565e+6;%5.9565e-6
a(8).c6=24.5325;
a(9).c6=0.6977;
a(10).c6=0.1108;
a(11).c6=0.015;
a(12).c6=403.43;


end






%deal_diff.m里面
function dy= deal_diff(t,y,a,istim,tstim)


dy=zeros(8,1);
P=[1.4,0.04,85,0.08,53,0.04,0.07,23,-0.04,0.1973,0.04,77,0.04,4,0.003,50,0.09,-82.3,-13.0287,-1e-7,0.07,1e-7]';
P1=P(1,1);
P2=P(2,1);
P3=P(3,1);
P4=P(4,1);
P5=P(5,1);
P6=P(6,1);
P7=P(7,1);
P8=P(8,1);
P9=P(9,1);
P10=P(10,1);
P11=P(11,1);
P12=P(12,1);
P13=P(13,1);
P14=P(14,1);
P15=P(15,1);
P16=P(16,1);
P17=P(17,1);
P18=P(18,1);
P19=P(19,1);
P20=P(20,1);
P21=P(21,1);
P22=P(22,1);

Cm=-1;



x1=y(1);
m=y(2);
h=y(3);
j=y(4);
d=y(5);
f=y(6);
ca=y(7);
vm=y(8);




dy(1)=(1-x1)*(a(1).c1*exp(a(1).c2*vm)+a(1).c3*(vm+a(1).c4))/(exp(a(1).c5*(vm+a(1).c4))+a(1).c6)...
-x1*(a(2).c1*exp(a(2).c2*vm)+a(2).c3*(vm+a(2).c4))/(exp(a(2).c5*(vm+a(2).c4))+a(2).c6);
dy(2)=(1-m)*(a(3).c1*exp(a(3).c2*vm)+a(3).c3*(vm+a(3).c4))/(exp(a(3).c5*(vm+a(3).c4))+a(3).c6)...
-m*(a(4).c1*exp(a(4).c2*vm)+a(4).c3*(vm+a(4).c4))/(exp(a(4).c5*(vm+a(4).c4))+a(4).c6);
dy(3)=(1-h)*(a(5).c1*exp(a(5).c2*vm)+a(5).c3*(vm+a(5).c4))/(exp(a(5).c5*(vm+a(5).c4))+a(5).c6)...
-h*(a(6).c1*exp(a(6).c2*vm)+a(6).c3*(vm+a(6).c4))/(exp(a(6).c5*(vm+a(6).c4))+a(6).c6);


dy(4)=(1-j)*(a(7).c1*exp(a(7).c2*vm)+a(7).c3*(vm+a(7).c4))/(exp(a(7).c5*(vm+a(7).c4))+a(7).c6)...
-j*(a(8).c1*exp(a(8).c2*vm)+a(8).c3*(vm+a(8).c4))/(exp(a(8).c5*(vm+a(8).c4))+a(8).c6) ;
dy(5)=(1-d)*(a(9).c1*exp(a(9).c2*vm)+a(9).c3*(vm+a(9).c4))/(exp(a(9).c5*(vm+a(9).c4))+a(9).c6)...
-d*(a(10).c1*exp(a(10).c2*vm)+a(10).c3*(vm+a(10).c4))/(exp(a(10).c5*(vm+a(10).c4))+a(10).c6);
dy(6)=(1-f)*(a(11).c1*exp(a(11).c2*vm)+a(11).c3*(vm+a(11).c4))/(exp(a(11).c5*(vm+a(11).c4))+a(11).c6)...
-f*(a(12).c1*exp(a(12).c2*vm)+a(12).c3*(vm+a(12).c4))/(exp(a(12).c5*(vm+a(12).c4))+a(12).c6);





ik1=P1*(exp(P2*(vm+P3))-1)/(exp(P4*(vm+P5))+exp(P6*(vm+P5)))+P7*(vm+P8)/(1-exp(P9*(vm+P8)));
ix1=x1*P10*(exp(P11*(vm+P12))-1)/exp(P13*vm);
ina=(P14*(m^3)*h*j+P15)*(vm-P16);
is=P17*d*f*(vm-P18-P19*log(ca));


dy(7)=P20*is+P21*(P22-y(7));
dy(8)=-(1/Cm)*(ik1+ix1+ina+is-istim) ;


end




怎么把deal_diff.m里面的中间变量ik1，ix1,ina,is传递出来？？
运行后提示：Warning: Failure at t=2.840902e-01. Unable to meet integration tolerances
without reducing the step size below the smallest value allowed
(1.009291e-15) at time t.
> In ode23t at 648
In second at 22