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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟 » beefly » 用matlab求动力学参数
51/1返回列表
查看: 1138  |  回复: 4
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

自制奶酪

铜虫 (小有名气)

[求助] 用matlab求动力学参数 已有2人参与

我不知道MATLAB怎么输入求动力学函数的方程，抄了木虫大佬的内容，请大家帮我把它弄到能运行出来，动力学参数拟合问题，求助大家~~~



function odes_fit
format long
clear all
clc
k0=[0 0 0 0 0 0 0 0];%参数初值
lb=[0 0 0 0 0 0 0 0];ub=[+inf +inf +inf +inf +inf +inf +inf +inf];%lb、ub:参数下限和上限
x0=[0 0 0 0 0];

data=...
    [0.5         1             1.5            2            3           4           5           6
          18.7    6.5    2.7    2.1   1.7   1.5   1.2   0.2;
          48.9    31.3   15.3   6.9   2.2   1.8   1.7   0.5;
          19.3    34.7   41.2   40.9   37.2   31.7   26.1   19.4;
          9.7    19.6    28.7   34.5   38.4   38.2   36.3   33.4;
          0.2    2.6    5.7    9.3   17.1   24.8   30.5   36.7;
];
x0=data(1,2:end);
tspan = [data(:,1)'];
yexp = [data(2:end,2) data(2:end,3) data(2:end,4) data(2:end,5) data(2:end,6) data(2:end,7) data(2:end,8)];
%使用函数fmincon()进行参数估计
[k,fval,flag]=fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.9f \n',k(1))
fprintf('\tk2 = %.9f \n',k(2))
fprintf('\tk3 = %.9f \n',k(3))
fprintf('\tk4 = %.9f \n',k(4))
fprintf('\tk5 = %.9f \n',k(5))
fprintf('\tk6 = %.9f \n',k(6))
fprintf('The sum of the squares is:%.le\n\n',fval)
k_fmincon=k;
%使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian]=...
    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值:\n'),output
%以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0=k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);
ci=nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值:\n')
output
figure(1)

ts=0(max(tspan)-min(tspan))/100):max(tspan);
[ts, ys] = ode45(@KineticsEqs,ts,x0,[],k);
yy = [data(:,2) data(:,3) data(:,4) data(:,5) data(:,6)];
figure(1)
plot(ts,ys(:,1),'b',tspan,yy(:,1),'bo');
figure(2)
plot(ts,ys(:,2),'r',tspan,yy(:,2),'ro');
figure(3)
plot(ts,ys(:,3),'k',tspan,yy(:,3),'ko');
figure(4)
plot(ts,ys(:,4),'g',tspan,yy(:,4),'ko');
figure(5)
plot(ts,ys(:,4),'g',tspan,yy(:,4),'ko');
figure(6)
plot(ts,ys(:,6),'m',tspan,yy(:,6),'mo');
%legend('C1的计算值','C1的实验值','C2的计算值','C2的实验值','C3的计算值','C3的实验值','C4的计算值','C4的实验值','C5的计算值','C5的实验值','C6的计算值','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);
Xsim5=Xsim(:,5);
Xsim6=Xsim(:,6);

ysim(:,1) = Xsim1(2:end);
ysim(:,2) = Xsim2(2:end);
ysim(:,3) = Xsim3(2:end);
ysim(:,4) = Xsim4(2:end);
ysim(:,5) = Xsim5(2:end);
ysim(:,6) = Xsim6(2:end);



f = [(ysim(:,1)-yexp(:,1)) (ysim(:,2)-yexp(:,2)) (ysim(:,3)-yexp(:,3)) (ysim(:,5)-yexp(:,5)) ...
    (ysim(:,6)-yexp(:,6))];

function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
C1=C(1);C2=C(2);C3=C(3);C4=C(4);C5=C(5);C6=C(6);
k1=k(1);k2=k(2);k3=k(3);k4=k(4);k5=k(5);k6=k(6);k7=k(7);
dC1dt = -k1*C1-k7*C1;
dC2dt = k1*C1-k2*C2-k5*C2;
dC3dt = k2*C2-k3*C3-k6*C3;
dC4dt = k3*C3-k4*C5;
dC5dt = k4*C5;
dC6dt = k5*C2+k6*C3+k7*C1;

dCdt = [dC1dt;dC2dt;dC3dt;dC4dt;dC5dt;dC6dt]; @beefly
回复此楼

» 猜你喜欢
简单快乐
1楼 2021-01-05 20:42:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

独孤神宇

版主 (知名作家)

【答案】应助回帖

感谢参与，应助指数 +1
这个代码不完整

lsqnonlin 函数 ObjFunc4LNL 缺失
赞一下 回复此楼
数值计算
2楼2021-01-05 21:03:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

自制奶酪

铜虫 (小有名气)
引用回帖:
2楼: Originally posted by 独孤神宇 at 2021-01-05 21:03:26
这个代码不完整

lsqnonlin 函数 ObjFunc4LNL 缺失

您好，这个要怎么修改啊，求助
赞一下 回复此楼
简单快乐
3楼2021-01-05 21:07:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
自制奶酪: 金币+20, ★有帮助 2021-01-05 21:39:36
独孤神宇: 金币+5, 鼓励交流 2021-01-06 17:08:32
参考下1stOpt的计算结果：
Weighted Root of Mean Square Error (RMSE): 3.52215359203778
Weighted Sum of Squared Residual: 434.194807406662
Correlation Coef. (R): 0.967264353331916
R-Square: 0.935600329226609
Determination Coef. (DC): 0.698480780131114
F-Statistic: -0.866225958773538

Parameter                  Best Estimate
--------------------        -------------
k1        1.18206935637483
k2        1.23804230026633
k3        0.392064236557302
k4        0.237399839174299
k5        2.4090372987211E-18
k6        1.44475738382174E-24
k7        0.719127412563461
f Initial Value         19.2237430015159
赞一下 回复此楼
4楼2021-01-05 21:34:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

流浪在地球

新虫 (初入文坛)
您好

发自小木虫Android客户端
回复此楼
5楼2022-11-17 13:58:31
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 自制奶酪 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭