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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 呼叫版主，在线紧急求助，关于matlab中微分方程组参数拟合得问题！
北京石油化工学院2026年研究生招生接收调剂公告
51/1返回列表
查看: 1757  |  回复: 12
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

lf0612

木虫 (小有名气)

[求助] 呼叫版主，在线紧急求助，关于matlab中微分方程组参数拟合得问题！ 已有2人参与

参考版主“月只蓝”以前帖子中的代码，写了如下关于微分方程组中参数拟合得代码，进行了修改，结果有问题，因为是新手，代码的含义不全懂，加之时间急着要，所以特来求助，往各位高手多多帮助，在此先谢谢了！
  
function k1k2k3
format long
clear all
clc
tspan = [0 1 2 3 4 5 6 7 8 9 10 ];
x0=[1 1 1 1 ];
k0 = [1 1 1 1 1 1 1 1 1 1 1 ]*1e-4;   
lb = [0 0 0 0 0 0 0 0 0 0 0];
ub = [1 1 1 1 1 1 1 1 1 1 1 ]*1e5;


data=[0.012,0.073,0.089,0.092,0.316,0.498,1.810,2.314,0.679,0.063,0.044;
    0.005,0.084,0.093,0.270,0.406,0.642,0.814,1.026,0.422,0.053,0.014;
    0.037,0.638,0.969,2.057,20.419,39.655,69.107,34.099,5.657,0.697,0.346;
   0.008,0.224,0.234,0.347,1.860,2.854,3.786,1.883,0.601,0.157,0.079;
   ]';
size(data);
yexp = data;

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))
fprintf('\tk5 = %.16f \n',k(6))
fprintf('\tk5 = %.16f \n',k(7))
fprintf('\tk5 = %.16f \n',k(8))
fprintf('\tk5 = %.16f \n',k(9))

ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,x0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,x0,[],k);

figure(1),
plot(ts,ys(:,1),'b',tspan(2:end),yexp(:,1),'or'),legend('计算值','实验值','Location','best');

figure(2),
plot(ts,ys(:,2),'b',tspan(2:end),yexp(:,2),'or'),legend('计算值','实验值','Location','best');

figure(3),
plot(ts,ys(:,3),'b',tspan(2:end),yexp(:,3),'or'),legend('计算值','实验值','Location','best');

figure(4),
plot(ts,ys(:,4),'b',tspan(2:end),yexp(:,4),'or'),legend('计算值','实验值','Location','best');



function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
Xsim4=Xsim(:,4);
ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
size(ysim(:,1));
size(ysim(:,2));
size(ysim(:,3));
size(ysim(:,4));
size(yexp(:,1));
size(yexp(:,2));
size(yexp(:,3));
size(yexp(:,4));
f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3))  (ysim(:,4)-yexp(:,4)) ];




function dEdt = KineticsEqs(t,C,k)   
% ODE模型方程 %g(1)=k(4);g(2)=k(5);g(3)=k(6);m(1)=k(7);m(2)=k(8);m(3)=k(9);
CA=C(1);CB=C(2);CC=C(3);CD=C(4);
%dCAdt = k(1)*CA*(CB/2-1)-k(2)*CB+k(3)*CC;
%dCBdt=-g(1)*CA+g(2)*CB.*(1-CB/5)+g(3)*CC.*(1-CC/10);
%dCCdt=m(1)*CA.*(CA/2-1)-m(2)*CB-m(3)*CC;
dCAdt= -k(1)*CA-k(2)*CA;
dCBdt=k(1)*CA-k(3)*CB-k(4)*CB;
dCCdt=k(3)*CB-k(5)*CC;
dCDdt=k(2)*CA+k(4)*CB+k(5)*CC;
dEdt=[dCAdt;dCBdt;dCCdt;dCDdt];

在此感谢！
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2014-06-24 21:08:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lf0612

木虫 (小有名气)
引用回帖:
3楼: Originally posted by 月只蓝 at 2014-06-24 21:38:38
如果你一定要用MATLAB的话，我可以帮你把程序调通。
不过对于常微分方程拟合，建议你找一下dingd专家，他有高版本的1stopt软件，计算高效强大。也许MATLAB编了半天程序，拟合结果没有1stopt花十分钟计算的结果好。

谢谢你了，帮忙把我把程序调通吧 ，我对参数的准确度要求不是那么的高，不要差的太多就好了！
赞一下 回复此楼
4楼2014-06-24 21:47:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 13 个回答

lf0612

木虫 (小有名气)
原始的微分方程组如下：
dS/dt=-k1*S-k2*S
dE/dt=k1*S-k3*E-k4*E
dI/dt=k3*E-k5*I
dR/dt=k2*S+-k4*E+k5*I
赞一下 回复此楼
2楼2014-06-24 21:17:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
lf0612: 金币+10, ★有帮助, 5 2014-07-04 16:29:05
如果你一定要用MATLAB的话，我可以帮你把程序调通。
不过对于常微分方程拟合，建议你找一下dingd专家，他有高版本的1stopt软件，计算高效强大。也许MATLAB编了半天程序，拟合结果没有1stopt花十分钟计算的结果好。
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
3楼2014-06-24 21:38:38
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

今晚比较忙，今天才着手做，不知道来得及不。程序已经调通
CODE:
function k1k2k3k4k5
format long
clear all
clc
tspan = [0 1 2 3 4 5 6 7 8 9 10];
x0=[0.012 0.005 0.037 0.008];
%k0 = [-0.5 -0.8 -0.8 -1.5 2];  
k0 = [-0.87 -0.9 48.11 -47.80 -0.142];
lb = -[1 1 1 1 1]*1e9;
ub = [1 1 1 1 1]*1e5;


data=[0.073,0.089,0.092,0.316,0.498,1.810,2.314,0.679,0.063,0.044;
    0.084,0.093,0.270,0.406,0.642,0.814,1.026,0.422,0.053,0.014;
    0.638,0.969,2.057,20.419,39.655,69.107,34.099,5.657,0.697,0.346;
   0.224,0.234,0.347,1.860,2.854,3.786,1.883,0.601,0.157,0.079;
   ]';
size(data);
yexp = data;

[k,resnorm,residual,exitflag,output,lambda,jacobian] =lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.16f \n',k(1))
fprintf('\tk2 = %.16f \n',k(2))
fprintf('\tk3 = %.16f \n',k(3))
fprintf('\tk4 = %.16f \n',k(4))
fprintf('\tk5 = %.16f \n',k(5))

ts=0:1:max(tspan);

[ts ys]=ode45(@KineticsEqs,ts,x0,[],k);
[ttt XXsim] = ode45(@KineticsEqs,tspan,x0,[],k);

figure(1),
plot(ts,ys(:,1),'b',tspan(2:end),yexp(:,1),'or'),legend('计算值','实验值','Location','best');

figure(2),
plot(ts,ys(:,2),'b',tspan(2:end),yexp(:,2),'or'),legend('计算值','实验值','Location','best');

figure(3),
plot(ts,ys(:,3),'b',tspan(2:end),yexp(:,3),'or'),legend('计算值','实验值','Location','best');

figure(4),
plot(ts,ys(:,4),'b',tspan(2:end),yexp(:,4),'or'),legend('计算值','实验值','Location','best');



function f = ObjFunc(k,tspan,x0,yexp)           % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
Xsim1=Xsim(:,1);
Xsim2=Xsim(:,2);
Xsim3=Xsim(:,3);
Xsim4=Xsim(:,4);
ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
size(ysim(:,1));
size(ysim(:,2));
size(ysim(:,3));
size(ysim(:,4));
size(yexp(:,1));
size(yexp(:,2));
size(yexp(:,3));
size(yexp(:,4));
f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3))  (ysim(:,4)-yexp(:,4)) ];




function dEdt = KineticsEqs(t,C,k)   
% ODE模型方程 %g(1)=k(4);g(2)=k(5);g(3)=k(6);m(1)=k(7);m(2)=k(8);m(3)=k(9);
CA=C(1);CB=C(2);CC=C(3);CD=C(4);
%dCAdt = k(1)*CA*(CB/2-1)-k(2)*CB+k(3)*CC;
%dCBdt=-g(1)*CA+g(2)*CB.*(1-CB/5)+g(3)*CC.*(1-CC/10);
%dCCdt=m(1)*CA.*(CA/2-1)-m(2)*CB-m(3)*CC;
dCAdt= -k(1)*CA-k(2)*CA;
dCBdt=k(1)*CA-k(3)*CB-k(4)*CB;
dCCdt=k(3)*CB-k(5)*CC;
dCDdt=k(2)*CA+k(4)*CB+k(5)*CC;
dEdt=[dCAdt;dCBdt;dCCdt;dCDdt];

参数比较多，你再调调初值k0吧。
多参数的拟合，还是用1stopt软件为好，在此推荐专家dingd。
赞一下(1人) 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
5楼2014-06-25 09:55:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 13 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 285求调剂 +4 AZMK 2026-04-02 5/250 2026-04-03 01:06 by 啵啵啵0119
[考研] 材料调剂 +8 懒羊羊轻置玉臀 2026-04-02 8/400 2026-04-02 22:03 by liu823948201
[考研] 交通运输考试264分求工科调剂 +4 jike777 2026-04-02 4/200 2026-04-02 21:53 by zllcz
[考研] 一志愿山东大学，085600，344 +7 魏子per 2026-04-02 8/400 2026-04-02 21:12 by 百灵童888
[考研] 315分 085602 求调剂 +10 26考研上岸版26 2026-04-02 10/500 2026-04-02 20:45 by dongzh2009
[考研] 318求调剂，计算材料方向 +10 吸喵有害笙命 2026-04-01 11/550 2026-04-02 16:29 by oooqiao
[考研] 调剂 +3 好好读书。 2026-04-01 6/300 2026-04-02 15:49 by liumengping
[考研] 08生物与医药专硕初试346找调剂 +6 dianeeee 2026-04-01 7/350 2026-04-02 08:23 by guoweigw
[考研] 085601材料工程找调剂 +20 oatmealR 2026-03-29 21/1050 2026-04-01 21:00 by lijunpoly
[考研] 379求调剂 +3 ?苦瓜不苦 2026-04-01 3/150 2026-04-01 20:09 by 暮云清寒
[考研] 材料专业调剂 +5 啦啦啦哭 2026-03-31 6/300 2026-04-01 16:48 by JourneyLucky
[考研] 考研材料工程351分调剂 +5 整个好的 2026-03-31 5/250 2026-04-01 09:36 by topgun2009
[考研] 333求调剂 +4 阿科逸 2026-03-31 4/200 2026-04-01 09:11 by jp9609
[考研] 085601 329分调剂 +6 yzsa12 2026-03-31 6/300 2026-03-31 15:23 by yanflower7133
[考研] 0703化学 +20 妮妮ninicgb 2026-03-27 20/1000 2026-03-31 13:33 by 无际的草原
[考研] 334求调剂 +7 Trying] 2026-03-31 7/350 2026-03-31 12:33 by 无际的草原
[考研] 283求调剂（080500） +14 A child 2026-03-27 14/700 2026-03-30 12:06 by 探123
[考研] 356求调剂 +3 gysy?s?a 2026-03-28 3/150 2026-03-29 00:33 by 544594351
[考研] 298调剂 +3 jiyingjie123 2026-03-27 3/150 2026-03-27 11:57 by wxiongid
[考研] 一志愿吉大071010，316分求调剂 +3 xgbiknn 2026-03-27 3/150 2026-03-27 10:36 by guoweigw
信息提示
关闭
请填处理意见
关闭