решение уравнения шредингера на языке pascal

[Den9713]

люди! помогите пожалуйста решить уравнение шредингера на языке pascal. Очень надо

[Circiter]

Cf. метод Нумерова. Может быть и примеры реализаций найдете.

[Circiter]

Чтобы вам было легче написать программульку, давайте рассмотрим самый простой случай.

Пусть в одномерном пространстве частица массы $\mu$ находится в потенциальном поле $U(x)$, а её состояние описывается волновой функцией $\psi(x)$ и не меняется со временем. Приравняв известную "полную" энергию $E$ к сумме кинетической и "внешней" потенциальной энергий (гамильтониану) получим одномерное стационарное уравнение Шрёдингера:

$$U(x)\psi(x)-\frac{\hbar^2}{2\mu}\frac{d^2}{dx^2}\psi(x)=E\psi(x).\eqno(1)$$

Для его численного интегрирования можно ввести регулярную сетку $x_1,\ \ldots,\ x_n$ с шагом $x_{i+1}-x_i=\delta$, а вторую производную волновой функции аппроксимировать например как $\psi''(x)\approx\delta^{-2}[\psi(x-\delta)+2\psi(x)+\psi(x+\delta)].$

Вычисление выражения $(1)$ в точках $x_2,\ \ldots,\ x_{n-1}$ даст $n-2$ линейных уравнений. Эти уравнения можно дополнить ещё двумя граничными условиями $\psi(x_1)=\psi(x_n)=0$, т.е., предположить, что частица прыгает в некоторой потенциальной яме. Полученную систему из $n$ уравнений нужно решить относительно $\psi(x_i)$.

Легко проверить, что исходная задача сводится к отысканию ядра матрицы $A$, определенной следующим образом: диагональные элементы $A_{i\ i}$ равны $U(x_i)-E+\hbar^2/(\mu\delta^2)$, элементы $A_{i\ i-1}$ и $A_{i\ i+1}$, т.е., элементы, расположенные рядом с диагональю равны $-\hbar^2/(2\mu\delta^2)$, а остальные элементы забиты нулями.

Вот, теперь можно решить такую систему, i.e., $A\psi=0$, записанную в матричной форме, и получить искомую собственную функцию гамильтониана (de facto, набор значений-отсчетов $\psi(x_i)$). Но это для очень простого случая, на самом деле все несколько сложнее.

Вообще, думаю, вам стоит уточнить условие задачи (ну напишите хотя-бы уравнение, которое вам нужно решить; укажите вид внешнего поля, граничные условия, etc), тогда может-быть кто-нибудь чего-нибудь да и подскажет дельного... Кроме того, неплохо было бы привести пример потенциала, для которого известно точное аналитическое решение; оно пригодится для тестирования вашей программки.

P.S.: Я не физик, почти ничего в мат. физике не понимаю, так-что если чепухи нагородил, не обращайте внимания. Позже, попробую, если это будет по-прежнему актуально, набросать таки что-нибудь на паскале, а пока лень...

[Circiter]

Придумал ещё одно сумасшедшее упрощение. На этот раз будут не только формулы, но и кусочки кода.

Выразим вторую производную волновой функции из уравнения $(1)$ предыдущего поста как $\psi''(x)=-2\mu\hbar^{-2}\left[E-U(x)\right]\psi(x)$ и введем для неё соответствующую переменную SecondDerivative. Коэффициент $-2\mu\hbar^{-2}$ обозначим как Mass, $E$ обозначим как Energy, а потенциалу $U(x)$ сопоставим функцию Potential(x). Кроме этого, нам понадобятся переменные FirstDerivative для $\psi'(x)$ и State для самой волновой функции $\psi(x)$. Переменная x хранит текущую точку пространства.

Теперь выражение для $\psi''(x)$ мы можем записать на паскале:

Используется синтаксис Pascal

    SecondDerivative:=Mass*(Energy-Potential(x))*State;
 

Производная $\psi'(x)$ может быть получена интегрированием $\psi''(x)$: $\psi'(x)=\int\psi''(x)dx$. Этот интеграл может быть приближен суммой $\psi'(x)\approx\sum_i\psi''(x_i)\delta$, где $\delta$ по-прежнему обозначает шаг сетки $\{x_i\}$ и в коде будет соответствовать константе dx (на каждой итерации происходит переход к слудующей точке по формуле x:=x+dx). Такое приближение сразу же можно перевести на паскаль:

Используется синтаксис Pascal

    FirstDerivative:=FirstDerivative+SecondDerivative*dx;
 

По аналогии заодно проинтегрируем ещё и эту производную, получив наконец-то саму $\psi(x)$:

Используется синтаксис Pascal

    State:=State+FirstDerivative*dx;
 

Реальную плотность вероятности обнаружения частицы в текущей точке можно получить просто возведя State в квадрат.

Все это крутим теперь в цикле, а цикл этот помещаем ещё в один цикл в котором перебираем подходящие значения переменной Energy (достаточно двигаться с единичным шагом). В общем, должно получиться что-то такое (псевдокод):

Используется синтаксис Pascal

const
    InitialPosition= ... ; {-10}
    dx=0.001;
    InitialEnergy= ... ; { 100 }
    InitialState= ... ; { -0.001 }
    MaximalState= ... ; {10000}
    InitialDerivative= ... ; { 0.001 }
    Mass= ... ; { 0.001 }

    ...

procedure Integrate(Energy: integer);
var
    x: real;
begin
    x:=InitialPosition;
    FirstDerivative=InitialDerivative;
    State=InitialState;

    repeat
        x:=x+dx;

        SecondDerivative=Mass*(Energy-Potential(x))*State;
        FirstDerivative=FirstDerivative+SecondDerivative*dx;
        State=State+FirstDerivative*dx;

        Yield(x, State*State);
    until State>MaximalState;
end.

procedure Solve;
var Energy: integer;
begin
    Energy:=InitialEnergy;

    repeat
        Integrate(Energy);
        Inc(Energy);
    until ... ;
end;
 

Здесь процедура Yield(x, p) может использоваться например для рисования очередной точки графика амплитуды вероятности. А основной процесс запускается процедурой Solve (программулька должна "подвиснуть" на "правильном" значении энергии). Осталось придумать, что должно быть на месте многоточий (в комментариях указаны возможные значения некоторых констант).

Функция Potential(x) может выглядеть примерно так (для частицы в настоящей потенциальной яме конечной глубины):

Используется синтаксис Pascal

function Potential(x: real): real;
const
    MaximalPotential= ... ; { Some large number, e.g., 10000}
    MinimalPotential= ... ; {1000}
    LeftBoundary= ... ; { Coordinate of the left wall, e.g., 0 }
    RightBoundary= ... ; { Right wall, e.g., 10 }
begin
    PotentialL=MaximalPotential;
    if (x>LeftBoundary)and(x<RightBoundary) then
        Potential:=MinimalPotential;
end;
 

P.S.: Не стоит эту писанину воспринимать шибко всерьез, просто используйте как пищу для размышлений или, на худой конец, как шаблон. Удачи