24小时热门版块排行榜    

Znn3bq.jpeg
查看: 3463  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[文学芳草园] 献血感触 +3 呀呀好傻 2026-05-19 3/150 2026-05-19 14:04 by jasswind
[有机交流] 反应很差,大量原料没有反应 5+3 Mr.Zot 2026-05-19 3/150 2026-05-19 13:56 by xtlilibin
[考博] 找博士生导师 +6 小代想上岸 2026-05-15 7/350 2026-05-19 10:22 by free_fisher
[教师之家] 上海大学实验技术岗位非升即走 +10 嘻嘻哈哈乐呵呵 2026-05-15 11/550 2026-05-19 10:03 by 嘻嘻哈哈乐呵呵
[基金申请] 评审有感 +10 popular289 2026-05-18 15/750 2026-05-19 09:18 by jurkat.1640
[考博] 26/27申博自荐-锂/钠电池方向 5+3 狗头军师. 2026-05-15 4/200 2026-05-19 09:10 by moonboat
[基金申请] 同样的基金本子,换个专家直接从C变A! +3 国自然国社科中 2026-05-19 3/150 2026-05-19 08:50 by Equinoxhua
[基金申请] 别被青基扩招骗了!26年科研内卷才刚刚开始 +3 国自然国社科中 2026-05-14 4/200 2026-05-19 08:48 by archvillain
[基金申请] 面上本子正文33页,违规吗?会被低分嘛? +8 1234567wang 2026-05-17 10/500 2026-05-18 18:52 by zzahkj
[基金申请] 国自然上会要求 +5 无名者登山 2026-05-18 9/450 2026-05-18 17:50 by BlakeReary
[基金申请] 今年审到国自然15份,谈谈感受 +16 国自然国社科中 2026-05-17 16/800 2026-05-18 14:58 by gy116024
[基金申请] 青C资助名额大幅增加! +12 西葫芦炒鸡蛋 2026-05-13 16/800 2026-05-18 10:02 by Equinoxhua
[文学芳草园] 半夜喝咖啡 +3 myrtle 2026-05-15 5/250 2026-05-18 01:03 by 小沈2018
[高分子] 本人最近太闲了,谁有问题可以提,每天会统一回复 +9 一切都是空工 2026-05-12 20/1000 2026-05-16 19:52 by Equinoxhua
[有机交流] 如何实现卤原子转化 +3 BT20230424 2026-05-15 5/250 2026-05-16 16:20 by czyzsu
[文学芳草园] 风把牡丹吹跑了 +5 myrtle 2026-05-12 9/450 2026-05-15 15:27 by myrtle
[教师之家] 教学课件你会给同学吗 +8 硕士研究生吗 2026-05-13 8/400 2026-05-14 22:23 by 常规沥青
[考博] 26应届毕业生考博求助 +3 wo一定上岸 2026-05-13 3/150 2026-05-14 21:47 by 明海天涯
[考博] 材料类只有一篇综述能申博么 +4 乐逍遥谷 2026-05-13 4/200 2026-05-14 12:05 by zhyzzh
[论文投稿] 求助大佬sci投稿哪个好中 +3 江沅188 2026-05-12 4/200 2026-05-13 14:35 by 江沅188
信息提示
请填处理意见