该微分方程组的解及绘图方法这是我自己编的程序,不知道为什么不能成功,求教>> syms s1 s2 x1 x2 t>> d_equt1='Ds1=(0.25*405+0.25*s2-0.44*s1)/(3+0.25*t)-5.654*10^(-6)*x1*s1';>> d_equt2='Ds2=(0.445*s1-0.25*s2/(9+0.25*t)-1.662
来源:学生作业帮助网 编辑:六六作业网 时间:2024/12/30 02:18:49
该微分方程组的解及绘图方法这是我自己编的程序,不知道为什么不能成功,求教>> syms s1 s2 x1 x2 t>> d_equt1='Ds1=(0.25*405+0.25*s2-0.44*s1)/(3+0.25*t)-5.654*10^(-6)*x1*s1';>> d_equt2='Ds2=(0.445*s1-0.25*s2/(9+0.25*t)-1.662
该微分方程组的解及绘图方法
这是我自己编的程序,不知道为什么不能成功,求教
>> syms s1 s2 x1 x2 t
>> d_equt1='Ds1=(0.25*405+0.25*s2-0.44*s1)/(3+0.25*t)-5.654*10^(-6)*x1*s1';
>> d_equt2='Ds2=(0.445*s1-0.25*s2/(9+0.25*t)-1.6625*10^(-5)*x2*s2*(405-s1)/(201+405-s1)';
>> d_equt3='Dx1=(0.25*x2-0.44*x1-0.25*x1)/(3+0.25*t)';
>> d_equt4='Dx2=(0.44*x1-0.25*x2-0.25*x2)/(9+0.25*t)';
>> condit1='s1(0)=30';
>> condit2='s2(0)=30';
>> condit3='x1(0)=3000';
>> condit4='x2(0)=3000';
>> [S1,S2,X1,X2,'t']=dsolve(d_eq1,d_eq2,d_eq3,d_eq4,condit1,condit2,condit3,condit4,'t')
该微分方程组的解及绘图方法这是我自己编的程序,不知道为什么不能成功,求教>> syms s1 s2 x1 x2 t>> d_equt1='Ds1=(0.25*405+0.25*s2-0.44*s1)/(3+0.25*t)-5.654*10^(-6)*x1*s1';>> d_equt2='Ds2=(0.445*s1-0.25*s2/(9+0.25*t)-1.662
你输入错了:
(1)第三行的第一个括号是多余的;
(2)所有的d_equt1, ...名称都和最后一行的d_eq1, ...不符,必须更正过来;
另外,没必要syms这些变量,然后还在最后一行多写两个't',matlab默认的自变量名称就是t,不怕它不认识.
算出来是没解析解的,于是改用ode45来求数值解并画图.先写如下的m-file,保存为myfun.m到当前目录下:
function y = myfun(t,x)
s1=x(1);s2=x(2);x1=x(3);x2=x(4);
y1=(0.25*405+0.25*s2-0.44*s1)/(3+0.25*t)-5.654*10^(-6)*x1*s1;
y2=0.445*s1-0.25*s2/(9+0.25*t)-1.6625*10^(-5)*x2*s2*(405-s1)/(201+405-s1);
y3=(0.25*x2-0.44*x1-0.25*x1)/(3+0.25*t);
y4=(0.44*x1-0.25*x2-0.25*x2)/(9+0.25*t);
y=[y1;y2;y3;y4];
之后,在命令窗口输入:
[t,x]=ode45(@myfun,[0:0.01:20],[30 30 3000 3000]);
s1=x(:,1);s2=x(:,2);x1=x(:,3);x2=x(:,4);
plot(t,s1,t,s2,t,x1,t,x2)
hleg1 = legend('s1','s2','x1','x2');
就是解方程和画图了.解完的变量储存在t和x里,x的第一列对应s1,第二列对应s2,以此类推.最后一行是在写图例,方便你看哪根线是哪个函数.注意,第一行的[0:0.01:20]是我自己设定的画图区间和步长,你自己可以自行调整,0.01是图上每个点的间隔,越小图越精确;0和20是t的起点和终点.你要是想画相图(这里顶多画三维的),比如你画s1和s2的,就是
plot(s1,s2)
画s1,s2,x1的,就是
plot3(s1,s2,x1)