| 查看: 1526 | 回复: 4 | |||
[求助]
用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))]; |
回复此楼» 收录本帖的淘帖专辑推荐
 常微分方程 |
» 猜你喜欢
青基已中
已经有3人回复
-
散金币祈福
已经有186人回复
-
物理学I论文润色/翻译怎么收费?
已经有118人回复
-
求助，如何提取ELK的rt-TDDFT在某一时刻的自旋密度分布和电子密度分布
已经有1人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
深圳大学应用技术学院招聘凝聚态物理博士后
已经有110人回复
-
深圳大学2027级光电信息工程/物理学研究生推免
已经有66人回复
-
求助中医院临床中医中标的青年标书学习，谢谢
已经有4人回复
-
浙大物院表面物理课题组长期招收普博及直博学生
已经有257人回复
独孤神宇
版主 (知名作家)
- 应助: 490 (硕士)
- 贵宾: 0.008
- 金币: 31030.3
- 散金: 802
- 红花: 122
- 沙发: 1
- 帖子: 5603
- 在线: 856.8小时
- 虫号: 3522474
- 注册: 2014-11-06
- 性别: GG
- 专业: 机械动力学
- 管辖: 计算模拟
2楼2018-05-11 09:35:55
3楼2018-05-11 12:06:27
独孤神宇
版主 (知名作家)
- 应助: 490 (硕士)
- 贵宾: 0.008
- 金币: 31030.3
- 散金: 802
- 红花: 122
- 沙发: 1
- 帖子: 5603
- 在线: 856.8小时
- 虫号: 3522474
- 注册: 2014-11-06
- 性别: GG
- 专业: 机械动力学
- 管辖: 计算模拟
4楼2018-05-11 15:41:36
5楼2018-05-18 13:56:38
10