大师作文网请问怎么用matlab求解既有微分方程又有一般方程的混合方程组啊
来源：学生作业帮 编辑：大师作文网作业帮 分类：数学作业 时间：2024/11/20 19:12:28
请问怎么用matlab求解既有微分方程又有一般方程的混合方程组啊
dy(1)=a1*t+b1*y(1);
y(2)=a2+b2*(t-y(1));
dy(3)=-a3*y(2)*(y(3)-t)/y(1) 请注意第二个方程没有微分,其中ab都是常数,想得到y(1)-t,y(2)-t,y(3)-t的图形
dy(1)=a1*t+b1*y(1);
y(2)=a2+b2*(t-y(1));
dy(3)=-a3*y(2)*(y(3)-t)/y(1) 请注意第二个方程没有微分,其中ab都是常数,想得到y(1)-t,y(2)-t,y(3)-t的图形
这个可以不当成方程组来处理,因为可以先解出第一个方程,再解出第二个方程,再解出第三个方程.另,dy(1)这样的表达方式看不出是对哪个自变量做的微分,并且没有提供初始条件或边界条件.
再问： 非常谢谢您啊。dy(1)，dy(3)都是对t的微分，初始条件知道，t范围也知道的。我可以解出第一个方程，第二个就不知道怎么解了，因为第一个方程只能得到y(1)-t的数值解，而y(2)跟t和y(1)都有关系啊，而y(3)中y(1)、y(2)都是随t变而改变，所以更不知道怎么解啊。我是matlab新手，恳请赐教啊！如果您能把相干程序附上感激不尽
再答： 嗯，客气了。我觉得你可能不知道怎么把第一个方程得到的解作为函数代入第二第三个方程，这个我在后面会补充说一下。先说将这个直接当作方程组处理吧。先编写方程组函数，例如： function dy = myfun(t,y) dy = zeros(3,1); % a column vector dy(1) = 2*t+2*y(1);% 这里的a1，b1参数我是随意取的，下面的也是 y(2) = 3+1*(t-y(1)); dy(3) = 2*y(2)*(y(3)-t)/y(1); 将它保存为myfun.m； 然后运行[T,Y] = ode23(@myfun,[0 10],[0 0 0]); 调用函数并得到解。Y矩阵的每一行分别对应要求的Y(1)、Y(2）、Y(3)； 至于要得到图形，比如y(1)-t，运行plot(T,Y(:,1)-T(:,1))即可。 =================================== 至于先求解了第一个方程后，解接下来的方程，这是另一种做法。 第二个方程其实不用解，就是个矩阵运算；第三个方程实际上相当于求解y'(t) + f(t)y(t) = g(t)类型的常微分方程，这里的f(t),g(t)由上方解已知，将它们代入第三个方程处理的方法是编写函数 function dydt = myfunn(t,y,T,f,T,g) f = interp1(T,f,t); % 对已知的f矩阵根据新的求解过程进行插值 g = interp1(T,g,t); % dydt = -f.*y + g 调用的话可以用[T Y] = ode45(@(t,y) myode(t,y,T,f,T,g),Tspan,IC)，Tspan和IC是自己定义的时间跨度和初始条件。 具体需要你自己调试了。 显然第一种方法更简单；虽然第二种方法有时候运算更有效率。
再问： 第一种方法很好用哦，非常感谢你。我的方程是这样的dy(1)=(0.558*abs(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))+20.996*(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))^0.5+128.852*(200*exp(-t/0.22)-y(1))+19.598)/(y(1)+10); y(2)=1036+159.12*(abs(200*exp(-t/0.22)-y(1)))^0.5; dy(3)=-3.744*y(2)*(y(3)-353+40*exp(-10*t))/(y(1)+10); y(1)-t是对的，可是y(2)-t画图怎么是直线呢？条件[0 0.4],[10 3229 780]
再答： 嗯，我刚发现我设置的初始值有问题，而且ode算法确实不支持混合方程。这样改一下myfun： function dy = myfun(t,y) dy = zeros(2,1); % a column vector dy(1)=(0.558*abs(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))+20.996*(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))^0.5+128.852*(200*exp(-t/0.22)-y(1))+19.598)/(y(1)+10); dy(2) = -3.744*(1036+159.12*(abs(200*exp(-t/0.22)-y(1)))^0.5)*(y(2)-353+40*exp(-10*t))/(y(1)+10); 求解的话调用[T,Y]=ode23(@myfun,[0,0.4],[10 780]);这样解出了Y(1)和Y(3)，当然Y(3)是用Y(2)表示的。然后再用Y(:,3)=Y(:,2); Y(:,2)=1036+159.12*(abs(200*exp(-T./0.22)-Y(:,1))).^0.5; 得到完整的Y矩阵，接下来就可以作图了。这是有重复计算的做法，应该不是最优的，也许把Y(2)设置成其他形式可以避免重复运算。
再问： 非常谢谢您啊。dy(1)，dy(3)都是对t的微分，初始条件知道，t范围也知道的。我可以解出第一个方程，第二个就不知道怎么解了，因为第一个方程只能得到y(1)-t的数值解，而y(2)跟t和y(1)都有关系啊，而y(3)中y(1)、y(2)都是随t变而改变，所以更不知道怎么解啊。我是matlab新手，恳请赐教啊！如果您能把相干程序附上感激不尽
再答： 嗯，客气了。我觉得你可能不知道怎么把第一个方程得到的解作为函数代入第二第三个方程，这个我在后面会补充说一下。先说将这个直接当作方程组处理吧。先编写方程组函数，例如： function dy = myfun(t,y) dy = zeros(3,1); % a column vector dy(1) = 2*t+2*y(1);% 这里的a1，b1参数我是随意取的，下面的也是 y(2) = 3+1*(t-y(1)); dy(3) = 2*y(2)*(y(3)-t)/y(1); 将它保存为myfun.m； 然后运行[T,Y] = ode23(@myfun,[0 10],[0 0 0]); 调用函数并得到解。Y矩阵的每一行分别对应要求的Y(1)、Y(2）、Y(3)； 至于要得到图形，比如y(1)-t，运行plot(T,Y(:,1)-T(:,1))即可。 =================================== 至于先求解了第一个方程后，解接下来的方程，这是另一种做法。 第二个方程其实不用解，就是个矩阵运算；第三个方程实际上相当于求解y'(t) + f(t)y(t) = g(t)类型的常微分方程，这里的f(t),g(t)由上方解已知，将它们代入第三个方程处理的方法是编写函数 function dydt = myfunn(t,y,T,f,T,g) f = interp1(T,f,t); % 对已知的f矩阵根据新的求解过程进行插值 g = interp1(T,g,t); % dydt = -f.*y + g 调用的话可以用[T Y] = ode45(@(t,y) myode(t,y,T,f,T,g),Tspan,IC)，Tspan和IC是自己定义的时间跨度和初始条件。 具体需要你自己调试了。 显然第一种方法更简单；虽然第二种方法有时候运算更有效率。
再问： 第一种方法很好用哦，非常感谢你。我的方程是这样的dy(1)=(0.558*abs(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))+20.996*(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))^0.5+128.852*(200*exp(-t/0.22)-y(1))+19.598)/(y(1)+10); y(2)=1036+159.12*(abs(200*exp(-t/0.22)-y(1)))^0.5; dy(3)=-3.744*y(2)*(y(3)-353+40*exp(-10*t))/(y(1)+10); y(1)-t是对的，可是y(2)-t画图怎么是直线呢？条件[0 0.4],[10 3229 780]
再答： 嗯，我刚发现我设置的初始值有问题，而且ode算法确实不支持混合方程。这样改一下myfun： function dy = myfun(t,y) dy = zeros(2,1); % a column vector dy(1)=(0.558*abs(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))+20.996*(200*exp(-t/0.22)-y(1))*(200*exp(-t/0.22)-y(1))^0.5+128.852*(200*exp(-t/0.22)-y(1))+19.598)/(y(1)+10); dy(2) = -3.744*(1036+159.12*(abs(200*exp(-t/0.22)-y(1)))^0.5)*(y(2)-353+40*exp(-10*t))/(y(1)+10); 求解的话调用[T,Y]=ode23(@myfun,[0,0.4],[10 780]);这样解出了Y(1)和Y(3)，当然Y(3)是用Y(2)表示的。然后再用Y(:,3)=Y(:,2); Y(:,2)=1036+159.12*(abs(200*exp(-T./0.22)-Y(:,1))).^0.5; 得到完整的Y矩阵，接下来就可以作图了。这是有重复计算的做法，应该不是最优的，也许把Y(2)设置成其他形式可以避免重复运算。