Рассмотрим соотношения между силой и энергией (работой). Для систем, в которых энергия сохраняется, можно сформулировать уравнения движения, используя скалярные величины, представляющие кинетическую и потенциальную энергию для частицы в силовом поле. Для неконсервативных сил интеграл работы следует рассчитывать для индивидуальных траекторий частиц с использованием линейного интеграла. Здесь рассматриваются оба случая.

Начнем с консервативного случая. Возьмем двухцентровый электростатический потенциал, который действует на валентный электрон в двухатомной молекуле (например, молекула ). (Кстати, совершенно аналогичным образом решается задача для космического корабля, астероида, маленькой планеты, движущейся в поле двух почти неподвижных звезд.) Итак, два протона находятся на фиксированном расстоянии 2d по оси z, один расположен в точке (0, 0, –d), другой – в точке (0, 0, d). На электрон в точке (x, y, z) действует сила притяжения, так как его заряд противоположен по знаку. Удобно выразить расстояние электрона от любого протона через его координаты (x, y, z), используя координаты r1 и r2 соответственно:

• restart
• r1:=sqrt(x^2+y^2+(z-d)^2);
• r2:=sqrt(x^2+y^2+(z+d)^2);

Потенциальная электростатическая энергия в единицах СИ задается через заряды протонов Q и электрона q, диэлектрическую проницаемость ε[0] и два расстояния:

• f[0]:=q*Q/(4*Pi*epsilon[0]);
• Usi:=-f[0]*(1/r1+1/r2);

Для построения графиков, а также для расчетов удобнее измерять энергию в единицах, в которых диэлектрическая проницаемость, единичный заряд и масштаб характерных длин включены в одну константу. Характерным масштабом длины можно выбрать d (или 2d), т. е. равновесное расстояние ядер в молекуле. В этих единицах имеем:

• U:=unapply(subs(d=1,-1/r1-1/r2),x,y,z)

Существуют различные способы визуализации этой функции, но важно, чтобы 2D-график был информативен. Для этого две координаты должны быть фиксированы, причем расстояния измеряются в единицах, кратных половине расстояния между ядрами d/2. Берутся сечения вдоль оси x (x = 0) с фиксированным y (рассматривается как параметр y0) и отображаются как функция от z. Чтобы избежать масштабирования, которое может понадобиться из-за сингулярности в потенциале, на графике отображается только нужный для задачи участок координат:

• P1:=plot(U(0,0,z),z=-3..3,-5..0,color=blue):
• P2:=plot(U(0,0.5,z),z=-3..3,color=red):
• P3:=plot(U(0,1,z),z=-3..3,color=green):
• with(plots):
• display({P1,P2,P3});

Из этого рисунка сразу для трех значений y0 (0, 0.5 и 1) получаются z-компоненты силы для электрона, расположенного в точке (0, y0, z). (Сам рисунок не показан. Получите его самостоятельно.). Принцип суперпозиции позволяет рассматривать эти потенциалы как одномерные потенциалы при движении с фиксированными x0 и y0, и, таким образом, сила просто задается как производная по z со знаком минус:

• P1f:=plot(-diff(U(0,0,z),z),z=-3..3,-5..5,color=blue,
  discont=true):
• P2f:=plot(-diff(U(0,0.5,z),z),z=-3..3,color=red):
• P3f:=plot(-diff(U(0,1,z),z),z=-3..3,color=green):
• display({P1f,P2f,P3f});

Чтобы избежать появления вертикальной линии в месте сингулярности, в график при x0 = 0 добавлена опция (как и положено делать в Maple). Обратите внимание, что знак силы изменяется и что положение точек равновесия (пересечение с нулем) зависит от значения y0. При y0 = 0 поведение качественно отличается и появляется особенность.

В традиционном подходе, чтобы отобразить полное трехмерное поведение, следовало бы перейти в цилиндрические координаты и использовать осевую симметрию по оси z. Но мы предпочитаем декартовы координаты. Это сохраняет ограничение x0 = 0, а потенциал будет отображаться в плоскости (yz), для чего используем график поверхности и контурную диаграмму:

• plot3d(U(0,y,z),y=-3..3,z=-3..3,view=-3..0,style=patchcontour,
  shading=zhue,axes=boxed,numpoints=3000,
  orientation=[-30,15]);

График поверхности содержит нарисованные на разной высоте контуры. Его можно превратить в контурную диаграмму, выбрав правильную проекцию (выбор углов θ, ). Используем в качестве альтернативы команду contourplot из пакета plots. Чтобы избежать сингулярности, выберем значения x0 = 0.1:

• contourplot(U(0.1,y,z),z=-3..3,y=-3..3,axes=boxed,
  scaling=CONSTRAINED,grid=[40,40],contours=[-4,-3,-2,-1.5,-1],
  filled=true);

Сила, действующая на частицу, определяется как градиент со знаком минус от соответствующей функции. Загружаем пакет linalg и получаем:

• with(linalg):
• F:=grad(-U(x,y,z),[x,y,z]);

Обратите внимание, что нельзя просто вычислить F: = – grad (U), так как это будет означать умножение вектора на скаляр (а именно на –1). Чтобы иметь возможность вычислять, следует применить evalm (эта процедура используется для умножения матриц и векторов).

Строим график для подтверждения предыдущего результата для z-компонента силы:

• plot(subs(x=0,y=0.5,F[3]),z=-3..3,-5..5,color=red);

Можно показать, что поле силы – безвихревое:

• curl(F,[x,y,z]);

Пример показал, что линейный интеграл не зависит от выбранного пути. Теперь выберем два пути в плоскости (yz), для которых потенциал показан на контурной диаграмме. Рассчитаем работу по перемещению из [0, 0, 0] до [0, 1, 1] двумя разными путями: один – по прямой линии, а другой – по дуге окружности с радиусом вокруг точки [0, 0.5, 0.5].

Чтобы записать линейный интеграл для вычисления работы, которую создает сила, действующая на частицу, проходящую по траектории r(t), используем параметрическое представление пути. Если путь дифференцируем по его параметру (например, по времени), то линейный интеграл можно свести к определенному интегралу Римана:

.

W = int(F . dr) = int(F(r(t)) . dr/dt, t = t0 .. t1).

Для прямого пути имеем:

• r_in:=[0,0,0]; r_fin:=[0,1,1];
• rSL:=evalm(r_in+t*(r_fin-r_in));
• rSL[1], rSL[2], rSL[3];

Без вызова evalm это не работает, а если rSL – вектор, команда ниже не работает. Сделаем так:

• vSL:=map(diff,rSL,t);
• whattype(rSL);

В параметрическом представлении прямой линии конечные точки записываются для t = 0 и t = 1 соответственно.

Предупреждение

Чтобы определенное интегрирование имело смысл и силы вычислялись правильно, т. е. только (!) вдоль выбранной траектории, в записи для силы нужно заменить x, y, z на значения x(t), y(t), z(t). Это наиболее важный аспект вычисления линейного интеграла, а также это наиболее распространенный источник ошибок).

В замене для всех компонент следует записать op(F), а не F:

• subs(x=rSL[1],y=rSL[2],z=rSL[3],op(F));

Используем процедуру dotprod из пакета linalg для вычисления скалярного произведения F(r(t)) и dr/dt, а затем интегрируем по t:

• W1:=int(dotprod(subs(x=rSL[1],y=rSL[2],z=rSL[3],
  op(F)),vSL),t=0..1);

Можно проверить, соответствует ли этот результат разнице между потенциальными энергиями, оцененными в двух конечных точках соответственно:

• U(op(r_fin))-U(op(r_in));
• evalf(%,4); evalf(W1,4);

Чтобы увидеть, как появилось первое из двух чисел, посмотрите на вклад в разность потенциалов:

• evalf(U(op(r_fin)),4); evalf(U(op(r_in)),4);

Количество работы по перемещению электрона из точки [0, 0, 0] в точку [0, 1, 1] отрицательно, что означает, что в среднем сила действовала против движения.

Попробуем более простые траектории, для которых вектор силы направлен одинаково с вектором положения, например переход от [0, 0, 0] к [0, 0, 1–c], т. е. от точки неустойчивого равновесия к одному из протонов, например к тому, который расположен в z = 1. Также посмотрим, что происходит, когда путь проходит через сингулярность…

Потенциальная энергия в начальном положении r_in=[0,0,0] более «притягательна», чем в конечном r_fin=[0,1,1]. Таким образом, разница в потенциальной энергии между конечной и начальной точками траектории положительна. Чтобы перейти от [0, 0, 0] к [0, 1, 1], надо затратить энергию. Это согласуется с предыдущим выводом о том, что работа отрицательна: чтобы переместить электрон из [0, 0, 0] в [0, 1, 1], надо затратить энергию. Для этого электрон может обмениваться кинетической энергией.

Теперь повторим расчет для круговой траектории:

• rC:=evalm(sqrt(2)/2*[0,cos(t),sin(t)]+[0,1/2,1/2]);
• rC[1];
• plot([rC[2],rC[3],t=Pi/4..5/4*Pi],scaling=CONSTRAINED);
• map(simplify,subs(t=Pi/4,op(rC)));
• map(simplify,subs(t=5*Pi/4,op(rC)));

При выбранной параметризации траектория начинается в 5/4π и заканчивается в π/4. Однако путь проходит через сингулярность в [0, 0, 1], и поэтому интегрирование выполняется в противоположном направлении:

• map(simplify,subs(t=-3*Pi/4,op(rC)));
• t0:=-3*Pi/4: t1:=Pi/4:
• vC:=map(diff,rC,t);
• W2:=int(expand(dotprod(subs(x=rC[1],y=rC[2],z=rC[3],
  op(F)),vC)),t=t0..t1)

Maple трудно взять интеграл: подынтегральное выражение (ниже) нетривиально. Если исключить особенность, то несложно интегрировать численно:

• arg:=expand(dotprod(subs(x=rC[1],y=rC[2],z=rC[3],op(F)),vC)):
• evalf(Int(arg,t=t0..t1),4);

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

Теперь, когда понятно, что потенциальная энергия позволяет вычислить работу, совершаемую безвихревой силой над частицами, посмотрим, как получить информацию о силе из контурной диаграммы. Поскольку эквипотенциальные контуры строятся в фиксированной плоскости, в ней будут видны только компоненты этой силы в этой плоскости. Сила задается градиентом потенциала (не забудьте про знак минус). Направление вектора градиента, возникающего в любой точке, ортогонально эквипотенциальному контуру, проходящему через эту точку. Пакет plots в Maple содержит команду, которая позволяет нарисовать набор стрелок на дискретной сетке для любой заданной потенциальной функции в двух измерениях. Чтобы избежать сингулярности, смотрим на плоскость (yz), определенную для x0 = 0.5:

• PL1:=gradplot(-U(0.5,y,z),z=-3..3,y=-2..2,axes=boxed,
  color=red,grid=[16,24],arrows=thick):
• PL2:=contourplot(U(0.5,y,z),z=-3..3,y=-2..2,axes=boxed,
  scaling=CONSTRAINED,grid=[40,40],contours=[-2,-1.5,-1,-0.5],
  filled=true,coloring=[white,blue]):
• display(PL1,PL2);

Стрелки дают хорошее представление об изменении двух составляющих силы в плоскости (yz). Информация о z-компоненте силы не передается.

График иллюстрирует тот факт, что интеграл работы не зависит от траектории. Видно, что при прохождении разных траекторий (конечные точки – разные) есть разнонаправленные векторы сил.

Проанализируйте движение электрона по пути, в котором он сначала ускоряется к протону, расположенному в [0, 0, d] в верхней полуплоскости, а затем разворачивается в направлении конечной точки [0, 1, 1], работая против силы, создаваемой полем (получится иначе, чем в предыдущем расчете прямого или кругового путей).

Таким образом, очень важна безвихревая природа консервативных сил.

В данном случае потенциал есть сумма двух центральных потенциалов относительно положений ядер в точках [0, 0, d] и [0, 0, –d]. Центральные потенциалы важны из-за их сферической симметрии (для таких потенциалов U зависит только от расстояния r, а не от сферических полярных углов θ и ). Такая симметрия подразумевает сохранение момента импульса как в классической, так и в квантовой физике.