应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » matlab 求指点 动力学方程拟合过程中导数的获取
北京石油化工学院2026年研究生招生接收调剂公告
111/2返回列表12下一页
查看: 3596  |  回复: 10
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

shidoudou

新虫 (初入文坛)

[求助] matlab 求指点 动力学方程拟合过程中导数的获取

本人matlab不通，勉强做了一个程序能够达到预期目的，可是老师想要中间运算的导数值，也就是dxdt的值，不知道如何填补表达式实现，望各位高手出手相助，不胜感激。
主程序：
format short
global Umax Ks Kp Y1 Y2 Uo Ko D Ux Kx Ku Ki;
Ki=0.0001;Ku=0.0006;
k0=[0.2 0.5 0.1 10 10 0.01 0.01 0.1 0.01 0.01]; %参数的初始值
x0=[0.0514 0.250 0]; %菌体浓度、底物浓度和产物浓度的初始值
t1=[0 6 12 18 24 30 36 42 48 54 60 66]';%发酵时间的实验数据
tspan=[0 6 12 18 24 30 36 42 48 54 60 66 ]';%时间间隔
%yexp实验数据[菌体浓度x1、底物浓度x2和产物浓度x3]
yexp=[[0.0514 0.0711 0.0997 0.1887 0.3478 0.6634 0.8711 1.1023 1.2336 1.4239 1.5849 1.6739];[0.25 0.2447 0.2366 0.2258 0.2196 0.203 0.1872 0.1698 0.1453 0.1124 0.0855 0.0608];[0 0.0012 0.0021 0.0039 0.0088 0.0137 0.0183 0.0268 0.0361 0.0438 0.0557 0.0643]]';
lb=[0 0 0 0 0 0.0216 0.001 0 0.001 0.001];ub=[1 10 10 50 50 2 1 1 0.1 0.1];%参数的下、上限
[k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc4LNL103,k0,lb,ub,[],x0,yexp);%非线性拟合
y1=[yexp(:,1)];y2=[yexp(:,2)];y3=[yexp(:,3)]';
[t4plot,x4plot]=ode45(@kineticseqs103,[0 100],x0,[],k);
plot(t1,y1,'b*',t4plot,x4plot(:,1),'k-'),xlabel('T(h)'),ylabel('Cx/（g/L）');
figure
plot(t1,y2,'g*',t4plot,x4plot(:,2),'k-'),xlabel('T(h)'),ylabel('Cs/(mmol/L)');
figure
plot(t1,y3,'r*',t4plot,x4plot(:,3),'k-'),xlabel('T(h)'),ylabel('Cp/(mmol/L)');
disp(k)
子程序：
1）
function f=ObjFunc4LNL103(k,x0,yexp)
tspan=[0 6 12 18 24 30 36 42 48 54 60 66];
[t1,x]=ode45(@kineticseqs103,tspan,x0,[],k);
y(:,1)=x(:,1);y(:,2)=x(:,2);y(:,3)=x(:,3);
f1=y(: ,1)-yexp(: ,1);f2=y(: ,2)-yexp(: ,2);f3=y(: ,3)-yexp(: ,3);
f=[f1*1 f2*5 f3*10];
2）
function dxdt=kineticseqs103(t,x,k) %模型方程
global Umax Ks Kp Y1 Y2 Uo Ko D Ux Kx Ku Ki
Umax=k(1); Ks=k(2); Kp=k(3); Y1=k(4); Y2=k(5); Ux=k(6); Kx=k(7); D=k(8); Uo=k(9); Ko=k(10);
  dxdt1=Umax*x(2)*x(1)*(1-Kp*x(3))/(Ks+x(2));
   if dxdt1<0
      dxdt1=0;
   end
  dxdt2=-1/Y1*dxdt1-Ux*x(2)*x(1)/(Kx+x(2))*(Ki/(Ki+x(3)));
  dxdt3=D*(-dxdt2-1/Y2*dxdt1-Uo*x(2)*x(1)/(Ko+x(2))*(Ku/(Ku+x(3))));
   if dxdt3<0
      dxdt3=0;
  end
  dxdt=[dxdt1;dxdt2;dxdt3];
只要能显示对应时间t1,这些点的dxdt值就可以了。谢谢
回复此楼

» 收录本帖的淘帖专辑推荐

matlab

» 本帖已获得的红花（最新10朵）

shidoudou

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2012-02-09 14:46:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
dbb627(金币+1): 欢迎交流 2012-02-10 18:25:29
最好用非Matlab代码的文本（图片）方式把问题描述一下。
赞一下 回复此楼
2楼2012-02-09 17:16:45
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

shidoudou

新虫 (初入文坛)
送鲜花一朵
非常感谢您哈。
我的本来的目的就是，建立了 微生物发酵的动力学方程，是三个微分方程，然后利用实验所得的数据就是yexp那组数据，非线性拟合微分方程里面的参数即可，以上的程序可以实现这个了。可是，老师让我能够显示这三个微分方程在对应的时间点t1那几个点处的导数值，我就不知道怎么样加语句了。不知道我表达清楚了没有～
赞一下 回复此楼
3楼2012-02-09 20:30:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hhucy

实习版主 (著名写手)

DOE锅炉工

【答案】应助回帖

感谢参与，应助指数 +1
dbb627(金币+1): 欢迎交流 2012-02-10 18:25:37
添加作图的命令，就能直观显示
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

shidoudou
人生那么多不确定，你怕什么
4楼2012-02-09 20:40:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

shidoudou

新虫 (初入文坛)
送鲜花一朵
引用回帖:
楼: Originally posted by hhucy at 2012-02-09 20:40:21:
添加作图的命令，就能直观显示

感谢哈。
能不能详细点致命咧，呵呵，我对此还真不通呢，可是老师说想要数值看一看，还得交差哈
赞一下 回复此楼
5楼2012-02-09 20:53:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师

【答案】应助回帖

感谢参与，应助指数 +1
dxdt=kineticseqs103(t,x,k) %模型方程

这不就是浓度梯度啊？
赞一下 回复此楼
6楼2012-02-10 15:58:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师

【答案】应助回帖

dbb627:编辑内容 2012-02-10 18:24
dbb627(金币+2): 感谢应助 2012-02-10 18:25:08
CODE:
[t,x]=ode45(@kineticseqs103,t1,x0,[],k);
n = length(x);
dxdt = zeros(size(x));
for i = 1:n
    dxdt(i,:) = kineticseqs103(t,x(i,:),k);
end
dxdt

[ Last edited by dbb627 on 2012-2-10 at 18:24 ]
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

shidoudou
7楼2012-02-10 16:30:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

715211229

木虫 (正式写手)
matlab命令窗口输入cftool，这个是拟合工具箱，拟合好后，最后一项分析里面有求一二阶导倒数的按钮
赞一下 回复此楼

» 本帖已获得的红花（最新10朵）

shidoudou
我是蜗牛
8楼2012-02-11 00:49:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

shidoudou

新虫 (初入文坛)
送鲜花一朵
引用回帖:
楼: Originally posted by change0618 at 2012-02-10 16:30:07:
CODE:
[t,x]=ode45(@kineticseqs103,t1,x0,[],k);
n = length(x);
dxdt = zeros(size(x));
for i = 1:n
    dxdt(i,:) = kineticseqs103(t,x(i,:),k);
end
dxdt

[ Last edited by dbb627 on  ...

谢谢哈
回复此楼
9楼2012-02-12 09:08:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

shidoudou

新虫 (初入文坛)
送鲜花一朵
引用回帖:
楼: Originally posted by 715211229 at 2012-02-11 00:49:28:
matlab命令窗口输入cftool，这个是拟合工具箱，拟合好后，最后一项分析里面有求一二阶导倒数的按钮

谢谢哈
回复此楼
10楼2012-02-12 09:08:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 shidoudou 的主题更新
111/2返回列表12下一页
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 285求调剂 +7 AZMK 2026-04-04 9/450 2026-04-06 00:06 by 永字号
[考研] 一志愿南昌大学，085600，344分求调剂 +6 调剂上岸玘 2026-04-05 6/300 2026-04-05 22:11 by 789风
[考研] 322求调剂 +3 嗯哼哼恒 2026-04-05 3/150 2026-04-05 19:52 by nepu_uu
[考研] 一志愿 南京航空航天大学 ，080500材料科学与工程学硕 +10 @taotao 2026-03-30 10/500 2026-04-05 17:57 by jj987
[考研] 266分，一志愿电气工程，本科材料，求材料专业调剂 +8 哇呼哼呼哼 2026-04-02 9/450 2026-04-05 17:14 by lbsjt
[考研] 一志愿南航，数一英一学硕317求调剂！！ +5 Acaciad 2026-04-04 5/250 2026-04-05 12:31 by 搏击518
[考研] 生物工程求调剂 +6 喜欢还是不甘心 2026-04-05 6/300 2026-04-05 10:28 by 唐沐儿
[考研] 一志愿电子科技大学085600材料与化工 329分求调剂 +10 Naiko 2026-04-04 10/500 2026-04-05 09:40 by sam3303
[考研] 070300化学学硕311分求调剂 +10 梁富贵险中求 2026-04-04 12/600 2026-04-05 09:37 by guoweigw
[考研] 调剂 +8 熊二想上岸 2026-04-04 8/400 2026-04-05 05:27 by houyaoxu
[考研] 323求调剂 +8 李佳乐1 2026-04-04 8/400 2026-04-04 22:26 by hemengdong
[考研] 可跨专业调剂 +3 周的得地 2026-04-04 6/300 2026-04-04 22:21 by barlinike
[考研] 调剂0855-288 +5 x熊二a 2026-04-03 5/250 2026-04-04 00:19 by 猪会飞
[考研] 288求调剂 一志愿哈工大 材料与化工 +39 洛神哥哥 2026-03-31 41/2050 2026-04-03 21:51 by qlm5820
[考研] 301求调剂 +14 A_JiXing 2026-04-01 14/700 2026-04-03 18:31 by ls刘帅
[考研] 11408，284分，二战真诚求调剂 +4 12.27 2026-04-02 4/200 2026-04-03 14:14 by dxiaoxin
[考研] 材料化工340求调剂 +5 jhx777 2026-03-30 5/250 2026-04-02 12:45 by smileboy2006
[考研] 求调剂 +4 DADA怪 2026-03-31 4/200 2026-04-01 14:30 by ZXlzxl0425
[考研] 358求调剂 +3 王向阳花 2026-03-31 3/150 2026-04-01 09:56 by zzchen2000
[考研] 调剂 +4 GK72 2026-03-30 4/200 2026-03-30 20:32 by dick_runner
信息提示
关闭
请填处理意见
关闭