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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » matlab M文件 给规定参数设个边界
81/1返回列表
查看: 918  |  回复: 7
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

cheng1378653

新虫 (初入文坛)

[求助] 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 ]
回复此楼

» 猜你喜欢
1楼2013-04-01 12:50:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

恩斯特

金虫 (小有名气)

【答案】应助回帖

★
感谢参与，应助指数 +1
fegg7502: 金币+1, 鼓励交流 2013-04-02 08:58:35
在nlinfit后面的第4个选项可以不填的吗？
赞一下 回复此楼
2楼2013-04-01 16:19:45
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lgycjpcqu

金虫 (正式写手)
★
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
赞一下(1人) 回复此楼
3楼2013-04-01 19:30:57
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cheng1378653

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by 恩斯特 at 2013-04-01 16:19:45
在nlinfit后面的第4个选项可以不填的吗？

我一直都这么用的，我以为还是初学者
赞一下 回复此楼
4楼2013-04-02 20:58:53
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)
lgycjpcqu的MATLAB代码就是正解。
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
5楼2013-04-04 16:37:22
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cheng1378653

新虫 (初入文坛)
引用回帖:
5楼: Originally posted by 月只蓝 at 2013-04-04 16:37:22
lgycjpcqu的MATLAB代码就是正解。

你能告诉我lgycjpcqu是什么吗？我以前没有听过这个。。。
赞一下 回复此楼
6楼2013-04-05 12:21:33
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)
引用回帖:
6楼: Originally posted by cheng1378653 at 2013-04-05 12:21:33
你能告诉我lgycjpcqu是什么吗？我以前没有听过这个。。。...

三楼的同学的名字，他提供的代码就是解决你的问题的。
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
7楼2013-04-05 14:05:26
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★
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
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
8楼2013-04-05 14:52:26
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 cheng1378653 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭