方势垒定态波函数程序(Matlab)

                     

贡献者: addis

预备知识 方势垒,Matlab 画图

   这里给出计算方势垒定态波函数的 Matlab 程序

图
图 1:FSB.m 画图结果

代码 1:FSB.m
% m 质量,V0 势垒,E 能量, L 势垒半宽, odd 奇函数 true 偶函数 false
% A 中间区间系数, C 右边区间 cos 系数, D 右边区间 sin 系数

function [psi,A,C,D] = FSB(x, E,L,m,V0,odd)
k = sqrt(2*m*E);
psi = zeros(size(x));
m1 = (x < -L); m2 = (abs(x) <= L); m3 = (x > L);
if (E == V0), E = V0 + eps; end

if E >= V0
    b = sqrt(2*m*(E - V0));
    if ~odd
        A1 = 1/sqrt(pi*(cos(b*L)^2+b^2*sin(b*L)^2/k^2));
        C1 = A1*(cos(k*L)*cos(b*L) + b/k*sin(k*L)*sin(b*L));
        D1 = A1*(-b/k*cos(k*L)*sin(b*L) + sin(k*L)*cos(b*L));
        A = A1; C = C1; D = D1;
        psi(m1) = C1*cos(k*x(m1)) - D1*sin(k*x(m1));
        psi(m2) = A1*cos(b*x(m2));
        psi(m3) = C1*cos(k*x(m3)) + D1*sin(k*x(m3));
    else
        B2 = 1/sqrt(pi*(sin(b*L)^2+b^2*cos(b*L)^2/k^2));
        C2 = B2*(cos(k*L)*sin(b*L) - b/k*sin(k*L)*cos(b*L));
        D2 = B2*(b/k*cos(k*L)*cos(b*L) + sin(k*L) * sin(b*L));
        A = B2; C = C2; D = D2;
        psi(m1) =-C2*cos(k*x(m1)) + D2*sin(k*x(m1));
        psi(m2) = B2*sin(b*x(m2));
        psi(m3) = C2*cos(k*x(m3)) + D2*sin(k*x(m3));
    end
else
    ka = sqrt(2*m*(V0-E));
    if ~odd
        A1 = 1/sqrt(pi*(cosh(ka*L)^2+(ka/k)^2*sinh(k*L)^2));
        C1 = A1*(cosh(ka*L)*cos(k*L) - ka/k*sinh(ka*L)*sin(k*L));
        D1 = A1*(cosh(ka*L)*sin(k*L) + ka/k*sinh(ka*L)*cos(k*L));
        A = A1; C = C1; D = D1;
        psi(m1) = C1*cos(k*x(m1)) - D1*sin(k*x(m1));
        psi(m2) = A1*cosh(ka*x(m2));
        psi(m3) = C1*cos(k*x(m3)) + D1*sin(k*x(m3));
    else
        B2 = 1/sqrt(pi*(cosh(ka*L)^2+(ka/k)^2*sinh(k*L)^2));
        C2 = B2*(sinh(ka*L)*cos(k*L) - ka/k*cosh(ka*L)*sin(k*L));
        D2 = B2*(sinh(ka*L)*sin(k*L) + ka/k*cosh(ka*L)*cos(k*L));
        A = B2; C = C2; D = D2;
        psi(m1) =-C2*cos(k*x(m1)) + D2*sin(k*x(m1));
        psi(m2) = B2*sinh(ka*x(m2));
        psi(m3) = C2*cos(k*x(m3)) + D2*sin(k*x(m3));
    end
end
end

   画图

代码 2:FSB_demo.m
% === 设置参数 =====
m = 1; V0 = 1; E = 2; L = 2;
% =================
x = linspace(-5, 5, 1000);
psi_e = FSB(x, E, L, m, V0, 0);
psi_o = FSB(x, E, L, m, V0, 1);
figure; plot(x, real(psi_e)); hold on; plot(x, real(psi_o));
xlabel x; legend('\psi_{even}', '\psi_{odd}');
结果如图 1

1. 行波

代码 3:FSB2.m
% m 质量,V0 势垒,E 能量, L 势垒半宽, odd 奇函数 true 偶函数 false
% A,B 左区间系数, C,D 中区间系数, E,F 右区间
% pm = 1: 从左入射; pm = -1: 从右边入射
function [psi, A,B,C,D,E,F] = FSB2(x, E, L, m, V0, pm)
k = sqrt(2*m*E);
psi = zeros(size(x));
m1 = (x < -L); m2 = (abs(x) <= L); m3 = (x > L);
if (E == V0), E = V0 + eps; end
if E >= V0
    b = sqrt(2*m*(E-V0));
    N = ((k^2+b^2)*sin(2*b*L)-2i*k*b*cos(2*b*L))*exp(-1i*k*L)...
         /sqrt(2*pi)/((k^2-b^2)^2*sin(2*b*L)^2 + 4*k^2*b^2);
    A = 1/sqrt(2*pi);
    B = N*(k^2-b^2)*sin(2*b*L)*exp(-1i*k*L);
    C = N*1i*k*(b+k)*exp(-1i*b*L);
    D = N*1i*k*(b-k)*exp(1i*b*L);
    E = N*2i*k*b*exp(-1i*k*L);
    F = 0;
    if pm < 0
        tmp = A; A = F; F = tmp;
        tmp = B; B = E; E = tmp;
        tmp = C; C = D; D = tmp;
    end
    psi(m1) = A*exp(1i*k*x(m1)) + B*exp(-1i*k*x(m1));
    psi(m2) = C*exp(1i*b*x(m2)) + D*exp(-1i*b*x(m2));
    psi(m3) = E*exp(1i*k*x(m3)) + F*exp(-1i*k*x(m3));
else
    ka = sqrt(2*m*(V0-E));
    N = ((k^2-ka^2)*sinh(2*ka*L)-2i*k*ka*cosh(2*ka*L))*exp(-1i*k*L)...
         /sqrt(2*pi)/((k^2+ka^2)^2*sinh(2*ka*L)^2 + 4*k^2*ka^2);
    A = 1/sqrt(2*pi);
    B = N*(k^2+ka^2)*sinh(2*ka*L)*exp(-1i*k*L);
    C = N*k*(1i*ka-k)*exp(-ka*L);
    D = N*k*(1i*ka+k)*exp(ka*L);
    E = N*2i*k*ka*exp(-1i*k*L);
    F = 0;
    if pm < 0
        tmp = A; A = F; F = tmp;
        tmp = B; B = E; E = tmp;
        tmp = C; C = D; D = tmp;
    end
    psi(m1) = A*exp(1i*k*x(m1)) + B*exp(-1i*k*x(m1));
    psi(m2) = C*exp(ka*x(m2)) + D*exp(-ka*x(m2));
    psi(m3) = E*exp(1i*k*x(m3)) + F*exp(-1i*k*x(m3));
end
end

   画图

代码 4:FSB2_demo.m
% === 设置参数 =====
m = 1; V0 = 0.9; E = 1.1; L = 2;
% =================
x = linspace(-5, 5, 1000);
psi = FSB2(x, E, L, m, V0, 1);
figure; plot(x, real(psi)); hold on; plot(x, imag(psi));
xlabel x;


致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 20 元,我们一周就能脱离亏损, 并在接下来的一年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利