24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

tandz

木虫 (正式写手)

应助: 0 (幼儿园)
金币: 1780.2
散金: 38
红花: 1
帖子: 449
在线: 703.8小时
虫号: 636814
注册: 2008-10-26
专业: 理论和计算化学

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

结合能E与距离r 已知. 请教各位，如何写个script拟合morse potential E(r) = D*(1-exp(a*(r-rm))^2的参数 D, a, rm呢? 本人几乎没有编程水平，请大家帮忙，非常感谢！

回复此楼

» 收录本帖的淘帖专辑推荐

精华

» 猜你喜欢

反铁磁体中的磁性切换：两种不同的机制已成功可视化已经有0人回复
求标准粉末衍射卡号 ICDD 01-076-1802 已经有0人回复
物理学I论文润色/翻译怎么收费? 已经有165人回复
新西兰Robinson研究所招收全奖PhD 已经有0人回复
石墨烯转移--二氧化硅衬底石墨烯已经有0人回复
笼目材料中量子自旋液体基态的证据已经有0人回复
数学教学论硕士可以读数学物理博士吗？已经有0人回复
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生已经有4人回复
澳门大学应用物理及材料工程研究院潘晖教授课题组诚招博士后已经有6人回复

1楼 2013-04-03 10:04:03

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

专家经验: +1059
计算强帖: 8
应助: 1712 (讲师)
贵宾: 8.888
金币: 68121.7
散金: 1938
红花: 443
沙发: 4
帖子: 4373
在线: 3291.4小时
虫号: 1122189
注册: 2010-10-14
专业: 宇宙学
管辖: 计算模拟区

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
感谢参与，应助指数 +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