在大神指点下　带着满满的诚意　重新发贴

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))];