24小时热门版块排行榜

返回列表

ankcean96

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 592.1
帖子: 161
在线: 4.9小时
虫号: 3092112
注册: 2014-03-27
专业: 化工冶金

[求助] 用matlab求动力学参数如何估计ｙ的初始值

在大神指点下　带着满满的诚意　重新发贴:victory::victory::victory:
matlab 求解动力学常微分方程组　　６个参数，初值和有物理学意义的区间如下：　　
k0 = [0.5575  0.0072  1  0.5136  0.01 0.2251];       % 参数初值
lb = [0.0001 0.0001 0.0001 0.1 0.01 0.1];                % 参数下限
ub = [50 1 1 1 1 1]; % 参数上限
常微分方程有五个，初始值如下：　
x0 = [1.4*10^8 2.043*10^4 1.8001*10^3 1212 0]; %y初值
按照下面写好的程序（y2初始值为2.043*10^4，y3初始值为1.8001*10^3），结果拟合的很好,但是实际上，y2和y3初始值是未知的。我只知道y2范围［1.8*10^4 ，1.4*10^7］，y3范围［1.8*10^3，1.8*10^5］，那么程序要怎么修改呢？

－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
format long
clear all
clc
lb = [0.0001 0.0001 0.0001 0.1 0.01 0.1];                % 参数下限
ub = [50 1 1 1 1 1]; % 参数上限
k0 = [0.5575  0.0072  1  0.5136  0.01 0.2251];       % 参数初值

x0 = [1.4*10^8 2.043*10^4 1.8001*10^3 1212 0]; %y初值
ExpData=[
31 3756.7 3217.9
59 2386.2 1060.9
90 7775.6 3470.9
120 24860.9 9981.9
151 35434.7 23046.0
181 34310 30359.4
212 26126.3 25344.2
243 11909.6 13215.4
273 10165.4 12080.2
304 8761.2 12249.1
334 7959.1 15276.8
365 6087.9 12667.9
];
t0=ExpData(:,1);
yexp =10* ExpData(:,2);  % yexp: CDC数据

% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@MSS,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('\tk3 = %.4f\n',k(3))
fprintf('\tk4 = %.4f\n',k(4))
fprintf('\tk5 = %.4f\n',k(5))
fprintf('\tk6 = %.4f\n',k(6))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fmincon = k;

%% =======绘图显示结果=======
tspan=[0:365];
[t x] = ode45(@li,tspan,x0,[],k);
yy(:,1) = x(:,4);
tt1=[sum(yy (1:31,:),1)
sum(yy (32:59,:),1)
sum(yy (60:90,:),1)
sum(yy (91:120,:),1)
sum(yy (121:151,:),1)
sum(yy (152:181,:),1)
sum(yy (182:212,:),1)
sum(yy (213:243,:),1)
sum(yy (244:273,:),1)
sum(yy (274:304,:),1)
sum(yy (305:334,:),1)
sum(yy (335:365,:),1) ]

plot(t0, yexp(:,1),'r.-')
hold on
plot(t0, tt1 (:,1),'bo-')
legend('y0', ' y1')
xlabel('t ')
ylabel('num')

function f = MSS(k,x0,yexp)          % fmincon优化
tspan=[0:365];
[t x] = ode45(@li,tspan,x0,[],k);
yy(:,1) = x(:,4);
tt1=[sum(yy (1:31,:),1)
sum(yy (32:59,:),1)
sum(yy (60:90,:),1)
sum(yy (91:120,:),1)
sum(yy (121:151,:),1)
sum(yy (152:181,:),1)
sum(yy (182:212,:),1)
sum(yy (213:243,:),1)
sum(yy (244:273,:),1)
sum(yy (274:304,:),1)
sum(yy (305:334,:),1)
sum(yy (335:365,:),1) ];
f = sum((log2(1+ yexp(:,1))-log2(1+ tt1(:,1))).^2);

%微分方程组
function dydt=li(t,y,k)
u=3.9139*10^-5;
dydt=[(u*(y(1)+y(2)+y(3)+y(4)+y(5))-k(1)*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))+k(2)*y(5)-u*y(1))
(k(1)*(y(3)+y(4))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5))-k(3)*y(2)-u*y(2))
(k(3)*(1-k(5))*y(2)-k(4)*y(3)-u*y(3))
(k(3)*k(5)*y(2)-k(6)*y(4)-u*y(4))
(k(4)*y(3)+k(6)*y(4)-k(2)*y(5)-u*y(5))];

回复此楼