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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » 呼叫版主，在线紧急求助，关于matlab中微分方程组参数拟合得问题！
131/2返回列表12下一页
查看: 1631  |  回复: 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的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

月只蓝

主管区长 (职业作家)

【答案】应助回帖

今晚比较忙，今天才着手做，不知道来得及不。程序已经调通
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的回帖
普通回帖

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的回帖

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的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
公式数据对不对啊？效果很差的。
赞一下 回复此楼
6楼2014-06-25 15:49:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lf0612

木虫 (小有名气)
引用回帖:
6楼: Originally posted by dingd at 2014-06-25 15:49:35
公式数据对不对啊？效果很差的。

谢谢了，数据有点问题，不过已经来不及了，所以最后手工试了几个参数，找了个最好的带进去了！再次感谢！
赞一下 回复此楼
7楼2014-07-04 16:31:45
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lf0612

木虫 (小有名气)
引用回帖:
5楼: Originally posted by 月只蓝 at 2014-06-25 09:55:08
今晚比较忙，今天才着手做，不知道来得及不。程序已经调通

function k1k2k3k4k5
format long
clear all
clc
tspan = ;
x0=;
%k0 = ;  
k0 = ;
lb = -*1e9;
ub = *1e5;


data=';
size(data);
yexp ...

谢谢了，数据有点问题，不过已经来不及了，所以最后手工试了几个参数，找了个最好的带进去了！再次感谢！
赞一下 回复此楼
8楼2014-07-04 16:31:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lf0612

木虫 (小有名气)
上网不方便，所以现在才回复二位，不好意思，谢谢了！
赞一下 回复此楼
9楼2014-07-04 16:32:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

小牧的影子

铁虫 (小有名气)
不错不错，学习了
回复此楼
学习，学习是立足之本。
10楼2015-06-04 21:47:45
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 lf0612 的主题更新
131/2返回列表12下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭