Question

MATLAB求微分方程组 并作图

dx/dt=a*b*c(1-c/y)
dy/dt=d*m(1-n*x/y+p)
其中，a,b,c,d,m,n,p都为常数。
且，a,b,c,d,m都是正整数，0<=n和p<=1；
a=？时可以使两式相等？

Answer 1

由于方程比较复杂，解析解不能用初等函数表示

只是要获得图像的话，用数值计算的办法可能更方便

fun=@(t,y) 1.44*(10^9)*(1-y).*exp(-109170./(8.314*t));

[T,Y] = ode23t(fun,500:600,0);

DY=fun(T,Y);

AX=plotyy(T,Y,T,DY);

set(get(AX(1),'Xlabel'),'String','T');

set(get(AX(1),'Ylabel'),'String','\alpha');

set(get(AX(2),'Ylabel'),'String','d\alpha/dT');

T是自变量，Y是变量也就是alpha

DY是，Y对T的导数

得到得到结果如下图

从图像看出，在T=500时，alpha=0

大概在T=600时，d(alpha)/dT趋向于0，alpha趋向于1不变

但是由于是数值解，在T不断增加的时候，d(alpha)/dT会在0附近振荡

所以T的取值不宜取得太大，这里取500到600之间

而采用ode23t函数，减少振荡

实际上，当T增大到一定值的时候，d(alpha)/dT趋向于0，

函数趋向于保持恒定值，所以后续的振荡是不合理的

取一定的区间如上图，已经可以很好地得到函数的变化趋势了

Answer 2

你是要画图还是求a=？时可以使两式相等？另外你说两式相等是说x(t) = y(t)还是dx/dt = dy/dt?
dsolve可以解微分方程组，ode45可以数值地解微分方程组，你需要数值解还是解析解？

追问

想知道二个方程的关系，也就是随着时间二者的x,y,值的变化规律，我想知道当常数a取值为多少时，可以使，y的值不在随着时间增加。
（当然想画图了，那样看最直观）但是我高数很烂，对那个软件也是很陌生。

追答

好的知道了。我用dsolve发现没有解析解（当然matlab说没解析解也不见得就真的解不出解析解，不过大部分情况下确实是没解析解的），就只能用数值解了。最终的变化规律还是要取决于你的常数a,b,c,d,m的数值。我写了个m函数文件如下，你先要在matlab编辑器窗口（不是命令窗口）输入这些：
function dz=myfun1(t,z)
global a b c d m n p
x = z(1);y = z(2);
dz(1) = a*b*c*(1-c/y);
dz(2) = d*m*(1-n*x/y+p);
dz = [dz(1);dz(2)];
这是定义一个新的函数。输入完后存在当前目录下，文件名为myfun1.m。完了你在命令窗口输入（我这里a,b,c,d这些常数的数值都是随便给的，t1,t2是画图的t的区间，x0,y0是解微分方程必须的初始条件，是x,y的初始点）：
global a b c d m n p
a=1;b=2;c=3;d=4;m=5;n=0;p=1;
t1=-2;t2=2;x0=1;y0=1;
[t,z]=ode45(@myfun1,[t1,t2],[x0,y0]);
subplot(2,1,1);
plot(t,z(:,1),t,z(:,2));
title('x,y关于t的函数图象');
subplot(2,1,2);
plot(z(:,1),z(:,2));
title('x,y的相图');
就完了。会画出两个图，第一个就是x,y关于t的图象，第二个是y关于x的图象，称为相图。变化规律你就看图说吧。
最后说你的问题，你想知道a取多少x,y不会随t变化。这就是要说两个导数始终为零，注意是始终为零。那么这就要求y = c 以及 x = (1+p)c/n恒成立，也就是都是常数函数才能达到你的要求，这个解和a没关系。我想这个应该不是你需要的答案，所以放最后说。上面的程序画出来的图象见我的图。需要审核一段时间才能看到的。
有问题继续问吧。我没上Hi的习惯，还是继续追问我吧。

追问

非常感谢你的细心指导，
如果限定那些参数，比如b=200;c=0.1;m=50;p=0.4;
可否找到a/d的与方程解的关系；就是这两个参数的比值，在某个范围时，可以使，y的值下降。或者不变。
（还有，就是怎么样知道图里那条线是X，哪个是y？）
这个是我改过常数值之后的图，想知道a/d对方程的影响。

追答

哪个x,哪个y，你改最后一部分就可以了。在
plot(t,z(:,1),t,z(:,2));
title('x,y关于t的函数图象');
中间插入一行：
legend('x','y')
即可。
你固定了b,c,m,p，但是n又不知道是多少，这种情况下要用程序解决你的问题还是得确定n才可能。原理就是尝试不同的n,a,d的数值的组合，解出在不同的数值组合下y的数值，然后画在一个坐标系里比较。具体你就要做循环了。比如尝试先固定n=2, a = 1, d = 1：
按上面的程序计算你会得到一组(x,y)，然后你输入：
y = z(:,2);
y2 = [];
y2 = [y2, y]
这样第一次计算得到的y的序列就存在y2里了。不要删除y2，删除其他变量，然后继续尝试
a = 1; d=2;
还是一样计算得到y
y = z(:,2);
y2 = [y2, y]
这样第二次计算得到的y的序列也存下来了。你做许多这样a,d的组合，算出矩阵y2后，然后作图
（我假设你做了四个a,d这样的组合）
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4))
legend('a=1, d=1','a=1,d=2','a=2,d=1','a=2,d=2')
画出来的图就是不同a,d组合下y随t变化的图形，你只有尝试不同的a,d的组合才可以，当然如果你是for a = 1:3, for d = 1:4，那就是12种组合，图形里有12根线，这个你再慢慢看吧。程序的解法就是很死且很繁琐的，不要对此寄予过高期望，加上你这里参数也很多，而且想要看两个变量对整个一条曲线的影响，这无疑是需要很多尝试的，计算也很费时。当然你这里dy/dt有特殊的形式，其实要找y不变或者下降就是找
dy/dt <= 0的条件，
这就是说d*m(1-n*x/y+p) <= 0，给定你这里dm都是正整数，那就是说要保证
1+p <= nx/y，所以你要看什么情况下y是否下降或者不变，就只要看左边这个不等式什么时候满足就可以了。虽然这个不等式里没有a和d，但是由于解出来的x,y是和a,d有关的，所以最终会和a,d有关系。这样解析地分析比数值计算的结果更加清晰，还不如这么做。