matlab 程序.rar_ode_常微分

function ptest(pval) % pval 是取值 16.5 17.5 17 % x0=0.1;y0=0.001;z0=0.001; %一组初始值 x0=1;y0=1;z0=2; %又一组初始值 q1=0.99;q2=0.99;q3=0.99; %分数导数的各个阶数（分别对应x,y,z变量） h=0.05;T=100; %h是步长，T为模拟时间 b=8/3;c=10; %微分方程中的系数 p=pval; %p初始化 %---- 下面三行是微分方程初始化 ----- x(1)=x0+h^q1*c*(y0-x0); y(1)=y0+h^q2*(x0*(p-z0)-y0); z(1)=z0+h^q3*(x0*y0-b*z0); t=h;n=1; %步数和时间 % ------ 开始迭代求解微分方程 -------- while t<=T %当还未超过模拟时间时，进行迭代 n=n+1; %步数计数 M1=0;M2=0;M3=0; %分数阶微分方程组各个函数离散化后的记忆项 % 下面的循环便是求解函数x，y，z各自的记忆项 for k=1:n-1 M1=M1+(-1)^(k)*gamma(q1+1)/(gamma(k+1)*gamma(q1-k+1))*(x(n-k)-x0); M2=M2+(-1)^(k)*gamma(q2+1)/(gamma(k+1)*gamma(q2-k+1))*(y(n-k)-y0); M3=M3+(-1)^(k)*gamma(q3+1)/(gamma(k+1)*gamma(q3-k+1))*(z(n-k)-z0); end t=t+h; % 迭代求解函数值, p值是确定的 x(n)=x0-M1+h^q1*c*(y(n-1)-x(n-1)); y(n)=y0-M2+h^q2*(x(n)*(p-z(n-1))-y(n-1)); z(n)=z0-M3+h^q3*(x(n)*y(n)-b*z(n-1)); end % 结束迭代求解过程 figure(1) %三维相轨迹 plot3(x,y,z);xlabel('x1'),ylabel('x2'),zlabel('x3'); title(sprintf('A0的三维相图(p1=%g)',pval)) figure(2) %x-y相轨迹 plot(x,y,'K');xlabel('x1'),ylabel('x2'); title(sprintf('A0的二维相图(p1=%g)',pval)) t=0:h:T; %时间离散 figure(3) %x(t)函数随时间变化图 subplot 221 plot(t,x,'K');xlabel('t'),ylabel('x1'); title(sprintf('A0的时域波形图(p1=%g)',pval)) subplot 222 %y(t)函数随时间变化图 plot(t,y,'K');xlabel('t'),ylabel('x2'); title(sprintf('A0的时域波形图(p1=%g)',pval)) subplot 223 %z(t)函数随时间变化图 plot(t,z,'K');xlabel('t'),ylabel('x3'); title(sprintf('A0的时域波形图(p1=%g)',pval))