[已解决] 根据模型建立微分方程组,并画出随时间变化的曲线

[复制链接]
要努力呀H 发表于 2021-1-20 17:06:25
本帖最后由 要努力呀H 于 2021-1-20 17:07 编辑

题目的要求是根据这个反应方案,画出和图A一样的图,但是我用我的程序编写出来之后和图A不一样,(其中一条是一样的)恳请各位指正,万分感谢。我的程序如下:
function sixpointfiveA()
%%本程序用于Mathematical Modelling Figure6.5A

x0=[0,0,0,0];
tspan=[0,1200];
[t,x]=ode23(@f,tspan,x0);
plot(t,x(:,1));
hold on;
plot(t,x(:,2),'rx');
hold on;
axis([0 1200 0 900]);
xticks(0:200:1200);
yticks(0:100:900);
title('Figure6.5A');
xlabel('Time(s)');
ylabel('Abundance(molecules per cell)');
function fy = f(t, x)
k1=0.002;k2=0.00001;k3=0.004;k5=0.01;k4=1;
if  100<=t&&t<=700
        L = 1;

    else
     L= 0;
end
fy = zeros(4, 1);
fy(1) = k2*(10000-x(1)-x(3)-x(4))*x(2)-k3*x(1);
fy(2) = k1*(4000-x(2))*L-k5*x(2);
fy(3) =-k4*x(4)*x(3)+k3*x(1);
fy(4)=-k4*x(4)*x(3)+k2*(10000-x(1)-x(3)-x(4))*x(2);
end
end

反应方案如图:
课本中的图案如图:

1.png
2.png

最佳答案


TouAkira 发表于 2021-1-21 01:14:32
要努力呀H 发表于 2021-1-20 12:44
您好,没有放出微分方程组确实是我的疏忽,但是关于k3取值的这个问题经过您的提醒,我发现这本书在印刷时 ...


啥??我就改了下k3的数值,得到的图看起来跟截图没什么差别。

  1. clear;
  2. x0=[0,0,0,0];
  3. tspan = [0,1200];
  4. [t,x]=ode23(@f,tspan,x0);
  5. plot( t, x( :, 1 ), 'k-', t, x( : , 2 ), 'r:', 'LineWidth', 2 );
  6. axis([0 1200 0 900]);
  7. xticks(0:200:1200);
  8. yticks(0:100:900);
  9. title('Figure6.5A');
  10. xlabel('Time(s)');
  11. ylabel('Abundance(molecules per cell)');
  12. function fy = f(t, x)
  13. k1=0.002;k2=0.00001;k3=0.11;k5=0.01;k4=1;
  14. if  100<=t&&t<=700
  15.         L = 1;
  16.     else
  17.      L= 0;
  18. end
  19. fy = zeros(4, 1);
  20. fy(1) = k2*(10000-x(1)-x(3)-x(4))*x(2)-k3*x(1);
  21. fy(2) = k1*(4000-x(2))*L-k5*x(2);
  22. fy(3) =-k4*x(4)*x(3)+k3*x(1);
  23. fy(4)=-k4*x(4)*x(3)+k2*(10000-x(1)-x(3)-x(4))*x(2);
  24. end
复制代码


[attach]318916[/attach]

5 条回复


TouAkira 发表于 2021-1-21 00:27:17
提问请把必要的信息给全了,你自己固然很清楚自己代码里面每个量对应图里哪个变量,别人不知道啊,直接在注释里面写明类似x( 1 )对应Ga、k1对应k_RL等信息,或者直接给出与你代码对应的微分方程组,对提问者来说不是什么做不到的事情吧,而且这种写注释的习惯,在你自己过几个月以后再看这份代码时,也能方便你快速准确地理解每一个环节。

光有一堆反应方程式,回答者还得再去推导一下各反应物的声称消耗微分方程。

这样强制要求回答者消耗本可以不必消耗的精力的时间去重复你已经做过一遍的工作,抬高了本来乐于帮你的人的时间成本,很大概率会导致原本有可能被解决的问题得不到解决——有人问路,大部分人都愿意花半分钟指一下方向,但如果让他先爬二十层楼再下来最后再花半分钟指路,换成是你,你还乐意吗,有这时间干点儿什么不好?

最后,我费了半天劲把各个物理量大致与微分方程对应上,然后发现,截图里面给的数值明明是 k_Gd0 = 0.11,  怎么到你写的自定义函数里面,就成了 k3=0.004 了?

这也忒粗心了吧?一个你自己完全可以检查出来的低级错误,非要用不直观的提问方式折腾掉回答者挺多时间,让我觉得特别憋屈,仿佛民警遇到个回不了家的小朋友,折腾半天才发现他没说清的家的地址就在眼前,而他所谓的回不了家是因为用错了钥匙开不了门。

但是来论坛提问的,普遍是高中以上学历了吧,数学和专业知识先不说,语文总归是经过九年义务教育锻炼的吧,总得学会清楚明确地表达问题吧。

再遇到问题,第一反应是先自己检查基本赋值、方程运算都写对没有;紧接着就是检查算法是否正确、语法是否符合版本要求;仍有问题再搜索论坛是否已有相似的被解决掉了的问题作为参考。上面这些都没能解决掉问题,才是给出必要信息并提问。

要努力呀H 发表于 2021-1-21 00:44:56
TouAkira 发表于 2021-1-21 00:27
提问请把必要的信息给全了,你自己固然很清楚自己代码里面每个量对应图里哪个变量,别人不知道啊,直接在注 ...

您好,没有放出微分方程组确实是我的疏忽,但是关于k3取值的这个问题经过您的提醒,我发现这本书在印刷时好像是出现了问题有两个图6.5,我在对照纸质版作图时k3的值的确为0.004(下面放出了那张图6.5),但无论如何,这确实是因为我的粗心所导致的不必要的问题,真的很抱歉。但是,刚刚我把值改为0.11时,图案仍和课本对不上,所以应该不是K3的问题,耽误了您的时间,我真的深感抱歉。
3.png

TouAkira 发表于 2021-1-21 01:14:32
要努力呀H 发表于 2021-1-20 12:44
您好,没有放出微分方程组确实是我的疏忽,但是关于k3取值的这个问题经过您的提醒,我发现这本书在印刷时 ...


啥??我就改了下k3的数值,得到的图看起来跟截图没什么差别。

  1. clear;
  2. x0=[0,0,0,0];
  3. tspan = [0,1200];
  4. [t,x]=ode23(@f,tspan,x0);
  5. plot( t, x( :, 1 ), 'k-', t, x( : , 2 ), 'r:', 'LineWidth', 2 );
  6. axis([0 1200 0 900]);
  7. xticks(0:200:1200);
  8. yticks(0:100:900);
  9. title('Figure6.5A');
  10. xlabel('Time(s)');
  11. ylabel('Abundance(molecules per cell)');
  12. function fy = f(t, x)
  13. k1=0.002;k2=0.00001;k3=0.11;k5=0.01;k4=1;
  14. if  100<=t&&t<=700
  15.         L = 1;
  16.     else
  17.      L= 0;
  18. end
  19. fy = zeros(4, 1);
  20. fy(1) = k2*(10000-x(1)-x(3)-x(4))*x(2)-k3*x(1);
  21. fy(2) = k1*(4000-x(2))*L-k5*x(2);
  22. fy(3) =-k4*x(4)*x(3)+k3*x(1);
  23. fy(4)=-k4*x(4)*x(3)+k2*(10000-x(1)-x(3)-x(4))*x(2);
  24. end
复制代码


ODEs65A.png

回复此楼

要努力呀H 发表于 2021-1-22 15:42:09
TouAkira 发表于 2021-1-21 01:14
啥??我就改了下k3的数值,得到的图看起来跟截图没什么差别。

很感谢您,我的MATLAB好像是反应有了一些问题,所以运行出来的图还是以前的,但是刚刚又验证了一遍,已经正确,再次表示感谢。

大大怪呀! 发表于 2021-1-26 08:32:30
可以使用龙格尔库塔法
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
热门教程
站长推荐
快速回复 返回顶部 返回列表