计算物理第1章-基本数学运算.zip_Lennard-Jones势_battle586

function semiclass clear clc close all format long tic EE=zeros(1,11); %for gamma=20:10:200; %对于特定分子的gamma值，如H2分子的gamma值为21.7，O2分子的gamma值为150 gamma=23; n=6; %振动能级数目; e0= -0.999; %简单搜索法启动点，LJ势函数（约化后）的最小值为-1 h1=0.0001; for k=1:n E(k)=qiugen(@(e)vh(e, k, gamma), e0, h1) %qiugen函数在下面子程序中定义,利用步长0.0001可以得到H2的6个束缚态,如果步长为0.001，则第六个束缚态无法求出，程序一直运行。 e0=E(k)+0.01*abs(E(k)); %改变初始猜测值（启动点）,用于求下一个振动能级 h1=0.01*abs(E(k)) [temp, jie(k,1), jie(k,2)] = vh(E(k),k, gamma); %获得r_in与r_out用于画图 end E2=E*4.747 %4.747eV为实验观测的位势深度 toc %EE((gamma-20)/10+1)=E(1); %EE1((gamma-20)/10+1)=E(2); %end %%对应于for gamma=循环 %plot(20:10:200,EE,'ro-',20:10:200,EE1,'bo-') figure;subplot(2,1,1); for k=1:n fplot(@(x)E(k), [jie(k,1),jie(k,2)]); %画振动能级 hold all end fplot(@(x)4*(x.^(-12)-x.^(-6)), [1 2]); %画势能函数曲线 hold off subplot(2,1,2) for k=1:n z=@(x)gamma*sqrt(abs(E(k)-4*(x.^(-12)-x.^(-6)))); %约化后的动量p，画相空间轨迹 fplot(@(x)[-z(x), z(x)], [jie(k,1) jie(k,2)]); hold on end %axis([1 2 -100 100]); function [f, r1,r2]=vh(e, n, gamma) ve=@(x)e-4*(x.^(-12)-x.^(-6)); %被积函数的平方，见教材（1.22）式 r1=qiugen(ve, 0.9, 0.001); %求r_in和r_out,0.9为简单搜索法启动点，0.001为步长 r2=qiugen(ve, 1.125, 0.001); %1.125启动点在势能面最低点处的自变量值2^(1/6)=1.1225附近 vee=@(x)(e-4*(x.^(-12)-x.^(-6))).^(1/2); %被积函数，见（1.22）式 f=(n-0.5)*pi-gamma*jifen(vee, r1, r2); %由于物理实际中n从0开始取值，此处从1开始取值，故为n-0.5 function jf=jifen(f, d, u, n) %f为被积函数，d,u分别为积分上下限，n为积分区间数 if nargin < 4 %nargin是用来判断输入变量个数的函数 n = 10000; %积分小区间数目 end h = (u-d)/n; %积分步长 t = d:h:u; ft = f(t); jf = (4*sum(ft(2:2:n))+2*sum(ft(3:2:n-1))+ft(1)+ft(n+1))*h/3; %复化simpson积分公式 function gen=qiugen(f, x0, h, delta) %求根函数，求解x_in和x_out两个点 if nargin < 4 %nargin是用来判断输入变量个数的函数 delta = 10^-7; %误差标准 end while h > delta %简单搜索法求根，可改成二分法，newton-raphson法等 if f(x0+h)*f(x0)>0 x0=x0+h; else h=h/2; end end % while h > delta %简单搜索结合二分法 % if f(x0+h)*f(x0)>0 % x0=x0+h; % else % a=x0; % b=x0+h; % while abs(a-b)>delta % x0=(a+b)/2; % if f(x0)*f(a)<0 % b=x0; % else % a=x0; % end % end % x0=(a+b)/2; % h=b-a; % end % end gen=x0;