| 查看: 1013 | 回复: 7 | ||||
[求助]
matlab M文件 给规定参数设个边界
|
|
Sample Text t=[0.0833 0.5 1 2 4 8]; c=[2.2505 1.6489 0.2789 0.2412 0.1077 0.1035]; 我用已知的两个函数模型拟合时出现了这种情况 （一）t=[0.0833 0.5 1 2 4 8]; c=[2.2505 1.6489 0.2789 0.2412 0.1077 0.1035]; myfun=inline('A(1)*exp(-A(2)*x)','A','x'); [A iter sse]=nlinfit(t,c,myfun,[1 1]) 结果：A = 2.6400 1.4099 iter = 6 sse = 0.2891 （二）：t=[0.0833 0.5 1 2 4 8]; c=[2.2505 1.6489 0.2789 0.2412 0.1077 0.1035]; myfun=inline('A(1)*exp(-A(2)*x)+A(3)*exp(-A(4)*x)','A','x'); [A iter sse]=nlinfit(t,c,myfun,[2.4 1.5 0.022 0.2]) 结果：A = 2.6266 1.4510 0.0227 -0.2002 iter = 19 sse = 0.2734 但是我的这个函数有个要求，就是参数A的值要大于零，这样的话A(4)就不符合了 （三）还有一种情况就是，用二的方式，但我给另外的随意的初值， 结果：A = 0.8571 1.4100 1.7829 1.4098 A(1)+A(3)=2.64，A(2)基本上等于A(4) iter = 6 sse = 0.2891 而这样的话就相当于把（一）的公式拆分成了两个一样模型的函数，与（一）没什么区别，我在用（二）计算式，我希望能够拟合出的结果中A(2)与A(4)不同，但是有还得是正值，所以我想给A(2)和A(4)设边界大于零 该怎么弄，我是matlab的菜鸟用户，请高手帮帮忙，谢谢！ [ Last edited by cheng1378653 on 2013-4-1 at 12:52 ] |
回复此楼» 猜你喜欢
帮我的英语口语老师找学生
已经有4人回复
-
计算三阶力常数用thirdorder扩胞时参数如何选择？临近原子数怎么取值
已经有7人回复
-
物理学I论文润色/翻译怎么收费?
已经有268人回复
-
2026年循环经济功能材料国际会议（ICFMCE 2026）
已经有0人回复
-
2026年第五届电气、电子与信息工程国际会议（ISEEIE 2026）
已经有0人回复
-
哈尔滨理工大学物理系招收物理学考研调剂
已经有5人回复
-
0702一志愿吉大B区求调剂
已经有5人回复
-
求调剂
已经有0人回复
-
0702一志愿吉大B区求调剂有论文
已经有0人回复
-
请问还有没有用Latex写文章的小伙伴们？
已经有0人回复
2楼2013-04-01 16:19:45
lgycjpcqu
金虫 (正式写手)
- 应助: 35 (小学生)
- 金币: 566.8
- 散金: 411
- 帖子: 371
- 在线: 176.9小时
- 虫号: 1327328
- 注册: 2011-06-20
- 性别: GG
- 专业: 电化学
★
fegg7502: 金币+1, 鼓励交流 2013-04-02 08:58:45
fegg7502: 金币+1, 鼓励交流 2013-04-02 08:58:45
|
function E=fconstrain(a,t,c) C=a(1).*exp(-a(2).*t)+a(3).*exp(-a(4).*t); E=c-C; 命令窗口 clear a0=[2.4 1.5 0.022 0]; LB=[-inf 0 -inf 0]; UB=[]; t=[0.0833 0.5 1 2 4 8]; c=[2.2505 1.6489 0.2789 0.2412 0.1077 0.1035]; options=optimset('lsqnonlin'); [A norm res ef]=lsqnonlin(@fconstrain,a0,LB,UB,options,t,c); cfit=A(1).*exp(-A(2).*t)+A(3).*exp(-A(4).*t); plot(t,c,'*',t,cfit,'rh') legend('原始数据','拟合数据') 结果是 A = 2.587176647067854 1.500899543303364 0.069042182557558 0.000000000000041  约束.jpg |
3楼2013-04-01 19:30:57
4楼2013-04-02 20:58:53
月只蓝
主管区长 (职业作家)
-
专家经验: +1059 - 计算强帖: 8
- 应助: 1712 (讲师)
- 贵宾: 8.888
- 金币: 68127.7
- 散金: 1938
- 红花: 443
- 沙发: 4
- 帖子: 4373
- 在线: 3291.4小时
- 虫号: 1122189
- 注册: 2010-10-14
- 专业: 宇宙学
- 管辖: 计算模拟区
5楼2013-04-04 16:37:22
6楼2013-04-05 12:21:33
月只蓝
主管区长 (职业作家)
-
专家经验: +1059 - 计算强帖: 8
- 应助: 1712 (讲师)
- 贵宾: 8.888
- 金币: 68127.7
- 散金: 1938
- 红花: 443
- 沙发: 4
- 帖子: 4373
- 在线: 3291.4小时
- 虫号: 1122189
- 注册: 2010-10-14
- 专业: 宇宙学
- 管辖: 计算模拟区
7楼2013-04-05 14:05:26
月只蓝
主管区长 (职业作家)
-
专家经验: +1059 - 计算强帖: 8
- 应助: 1712 (讲师)
- 贵宾: 8.888
- 金币: 68127.7
- 散金: 1938
- 红花: 443
- 沙发: 4
- 帖子: 4373
- 在线: 3291.4小时
- 虫号: 1122189
- 注册: 2010-10-14
- 专业: 宇宙学
- 管辖: 计算模拟区
【答案】应助回帖
★ ★ ★ ★
csgt0: 金币+2, 谢谢 2013-04-07 15:03:32
cheng1378653: 金币+2, ★★★很有帮助 2013-04-08 22:34:24
csgt0: 金币+2, 谢谢 2013-04-07 15:03:32
cheng1378653: 金币+2, ★★★很有帮助 2013-04-08 22:34:24
|
一下代码直接复制到m文件中即可。拟合结果以及相关性系数R^2见附图1。 %-------------------Start-------------------------------------------- function feixianxingnihe1 clear all;clc format long tspan=[0.0833 0.5 1 2 4 8]; %t的数据，在此输入 xexp=[2.2505 1.6489 0.2789 0.2412 0.1077 0.1035]; %c的数据，在此输入 k0=[2.6 1.45 0.022 0]; %这里设定A1~A4的初值 lb=-[1 0 1 0]*1e5; %分别设定A1~A4的取值下限，A2和A4已经设为大于等于0了 ub=[1 1 1]*1e5; %A1~A4的取值上限 %------------------------------------------------------------------------- % 使用函数lsqnonlin()进行参数估计 OPTIONS=optimset('MaxFunEvals',1000); [k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,k0,lb,ub,OPTIONS,tspan,xexp); ci = nlparci(k,residual,jacobian); %residual; fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\n\t参数 A1 = %.16f',k(1)) fprintf('\n\t参数 A2 = %.16f',k(2)) fprintf('\n\t参数 A3 = %.16f',k(3)) fprintf('\n\t参数 A4 = %.16f',k(4)) y=KineticsEqs(tspan,k); R2=1-sum((xexp-y).^2)./sum((xexp-mean(y)).^2); fprintf('\n\t相关系数之平方R^2 = %.16f',R2); figure plot(tspan,KineticsEqs(tspan,k),'b',tspan,xexp,'or'),legend('计算值','实验值','Location','Best') %------------------------------------------------------------------------- function f = ObjFunc(k,tspan,xexp) f=KineticsEqs(tspan,k)-xexp; %------------------------------------------------------------------------ function xt = KineticsEqs(t,k) xt=k(1)*exp(-k(2)*t)+k(3)*exp(-k(4)*t); %-----------------------------------The end ----------------------------- 在代码赋于初值的那一部分，可以根据各参数的物理意义输入具体数值，最后的拟合结果与初值关系很大。 另外用1stopt可以得到R^2=0.993409615207485的拟合结果： A1 -17577.3242768536 A2 99.2632248231063 A3 8.9625935335228 A4 3.38938238943609 以及若干组R^2=0.993409615207485的拟合结果： A1 8.96259352762256 A2 3.38938238846346 A3 -25330.6623331322 A4 103.649850688563 分析选定的方程形式可知，其实A1，A2，A3，A4具有对称性A(i)=8.96259352762256, 3.38938238846346是数值稳定的，剩余的两个A(i)只要一个数值较大，即会使得那一整项为0。 当然，R^2大的，不一定具有物理意义。  附图1.jpg |
8楼2013-04-05 14:52:26
