24小时热门版块排行榜    

查看: 3585  |  回复: 11

tandz

木虫 (正式写手)

[求助] 如何拟合morse potential的参数?

结合能E与距离r 已知. 请教各位,如何写个script拟合morse potential E(r) = D*(1-exp(a*(r-rm))^2的参数 D, a, rm呢? 本人几乎没有编程水平,请大家帮忙,非常感谢!
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
感谢参与,应助指数 +1
tandz: 金币+50, ★★★很有帮助, 太感谢了! 2013-04-04 17:03:15
tandz: 金币+20, ★★★★★最佳答案 2013-04-06 17:18:56
csgt0: 金币+2, 谢谢 2013-04-07 15:11:49
MATLAB代码可以吗?
假如已经获得了一组r,E的数据:
r:1 2 3 9
E:1 5 8 18
MATLAB程序如下,结果见附图1。
function feixianxingnihe
clear all;clc
format long

tspan=[1 2 3 9];    %r的数据,在此输入
xexp=[1 5 8 18];    %E的数据,在此输入

k0=[1 1 1];
lb=-[1 1 1]*1e5;
ub=[1 1 1]*1e5;


%-------------------------------------------------------------------------

% 使用函数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参数 D = %.16f',k(1))
fprintf('\n\t参数 a = %.16f',k(2))
fprintf('\n\t参数rm = %.16f',k(3))
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)*(1-exp(k(2).*(t-k(3)))).^2;

做方程拟合,1stopt代码简单易懂,1stopt代码如下:
Parameters D,a,rm;
Variable r,E;
Function E=D*(1-exp(a*(r-rm)))^2;
Data;
1 1
2 5
3 8
9 18
计算结果见附图2。   
可以看到两种方法计算结果相当接近。

附图1.jpg



附图2.jpg

MATLAB、MS小问题、普通问题请发帖求助!时间精力有限,恕不接受无偿私信求助。
5楼2013-04-04 16:27:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

dingd

铁杆木虫 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
tandz: 金币+10, ★★★很有帮助 2013-04-03 10:25:47
找个1stOpt,两分钟就能搞定。
2楼2013-04-03 10:15:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

tandz

木虫 (正式写手)

引用回帖:
2楼: Originally posted by dingd at 2013-04-03 10:15:16
找个1stOpt,两分钟就能搞定。

老板非要我写script, 没有办法啊
3楼2013-04-03 10:25:28
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

tandz

木虫 (正式写手)

自己顶,求助大家
4楼2013-04-03 19:21:49
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

tandz

木虫 (正式写手)

引用回帖:
5楼: Originally posted by 月只蓝 at 2013-04-04 16:27:04
MATLAB代码可以吗?
假如已经获得了一组r,E的数据:
r:1 2 3 9
E:1 5 8 18
MATLAB程序如下,结果见附图1。
function feixianxingnihe
clear all;clc
format long

tspan=;    %r的数据,在此输入
xexp ...

太感谢了!没学过MATLAB,我先试试。看起来效果非常好!
我学过一点点FORTRAN和C语言,不过现在都忘了,呵呵
6楼2013-04-04 17:06:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

tandz

木虫 (正式写手)

引用回帖:
5楼: Originally posted by 月只蓝 at 2013-04-04 16:27:04
MATLAB代码可以吗?
假如已经获得了一组r,E的数据:
r:1 2 3 9
E:1 5 8 18
MATLAB程序如下,结果见附图1。
function feixianxingnihe
clear all;clc
format long

tspan=;    %r的数据,在此输入
xexp ...

为什么我不能运行MATLAB程序呢,有几个错误
??? Error using ==> optim\private\lsqncommon
User supplied function failed with the following error:

Error using ==> feval
Undefined command/function 'ObjFunc'.

Error in ==> lsqnonlin at 147
[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...
7楼2013-04-06 16:06:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

引用回帖:
7楼: Originally posted by tandz at 2013-04-06 16:06:21
为什么我不能运行MATLAB程序呢,有几个错误
??? Error using ==> optim\private\lsqncommon
User supplied function failed with the following error:

Error using ==> feval
Undefined command/func ...

file-new-m文件或者function新建一个m文件,上述代码复制进m文件中,F5运行。我刚又测试过,可以运行。

» 本帖已获得的红花(最新10朵)

MATLAB、MS小问题、普通问题请发帖求助!时间精力有限,恕不接受无偿私信求助。
8楼2013-04-06 18:08:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

tandz

木虫 (正式写手)

送红花一朵
引用回帖:
8楼: Originally posted by 月只蓝 at 2013-04-06 18:08:39
file-new-m文件或者function新建一个m文件,上述代码复制进m文件中,F5运行。我刚又测试过,可以运行。...

运行了!谢谢帮助!
9楼2013-04-07 09:58:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

ytudou

银虫 (小有名气)

引用回帖:
5楼: Originally posted by 月只蓝 at 2013-04-04 16:27:04
MATLAB代码可以吗?
假如已经获得了一组r,E的数据:
r:1 2 3 9
E:1 5 8 18
MATLAB程序如下,结果见附图1。
function feixianxingnihe
clear all;clc
format long

tspan=;    %r的数据,在此输入
xexp ...

大神,你看他这个势函数形式没问题吗?morse势不是形如u=D{exp(-2α(r-ro))-2exp(-α(r-ro))}的样子吗?还是说我写出来的这个可以通过变形变成那个?请赐教
要给力啊,MS~LAMMPS~
10楼2013-10-11 14:37:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 tandz 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] fileCode有新解读? +7 Tide man 2026-08-08 13/650 2026-08-08 09:47 by zhengxianzb
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 HEQlVqMTIA7d 2026-08-07 5/250 2026-08-08 08:59 by 6vVgjDL4CnGu
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +3 HEQlVqMTIA7d 2026-08-07 4/200 2026-08-08 08:27 by 6vVgjDL4CnGu
[基金申请] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +3 KXLV3nuBVBY7 2026-08-07 4/200 2026-08-08 08:21 by 6vVgjDL4CnGu
[基金申请] 关于filecode +4 布布和一二 2026-08-07 7/350 2026-08-07 22:55 by zhanghaozhu
[基金申请] 娱乐 +6 Tide man 2026-08-03 6/300 2026-08-07 22:40 by 铁帽子农民
[基金申请] 化学口download_prp&fileCode的固定段好像这几天一直没变,有变的大神么? +3 Tide man 2026-08-07 4/200 2026-08-07 22:39 by Tide man
[基金申请] 关于代码变化问题,想知道的进来 +9 且听虎啸 2026-08-07 12/600 2026-08-07 21:53 by 理想三旬2025
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
[基金申请] 听说今天filecode变了 +24 布布和一二 2026-08-06 47/2350 2026-08-07 16:02 by zhiyanjiang
[基金申请] filecode变化情况 +6 布布和一二 2026-08-07 22/1100 2026-08-07 14:45 by 且听虎啸
[基金申请] 大家散了吧,后缀研究没有意义,别浪费时间了,过好目前的每一天,不要焦虑 +5 Tide man 2026-08-06 7/350 2026-08-07 13:11 by 医学老男孩
[基金申请] filecode +14 等待解的谜 2026-08-06 19/950 2026-08-07 12:20 by wlwhappy
[基金申请] filecode +8 布布和一二 2026-08-06 11/550 2026-08-06 20:41 by tangpu318
[基金申请] 8月时间戳变的,举个手。玩一下,释放压力 +9 archvillain 2026-08-04 11/550 2026-08-05 20:06 by wlwhappy
[基金申请] 好消息?这个有何含义??? +8 Tide man 2026-08-05 10/500 2026-08-05 16:14 by xmuxiaoyu
[考博] 【2027博士申请】纳米药物递送方向 20+3 13586093586 2026-08-03 4/200 2026-08-05 09:59 by lfy8008
[基金申请] 纯娱乐,不喜欢勿喷 +7 Tide man 2026-08-04 10/500 2026-08-04 15:10 by loufangrui
[基金申请] 面上提前没消息,有中的吗 +14 archvillain 2026-08-02 18/900 2026-08-04 14:42 by archvillain
[基金申请] 什么时候能放榜呀? +3 Jacob678 2026-08-03 3/150 2026-08-03 16:14 by gltch
信息提示
请填处理意见