Равновесное размещение одноименных зарядов на сфере
Арсений Слободюк, 2006г.

> with(linalg):with(plots):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

> assume(x1,real,y1,real,z1,real,x2,real,y2,real,z2,real,ftx,real,fty,real,ftz,real);

Составим систему уравнений, чтобы найти Ft - тангенциальную составляющую силы, действующей со стороны заряда

с координатами v1 на заряд с координатами v2. Заряды помещены на стержнях, закрепленных в начале координат.

> Ft:=vector(3,[ftx,fty,ftz]);

> v1:=vector(3,[x1,y1,z1]);

> v2:=vector(3,[x2,y2,z2]);

Сила Кулона. Для данной задачи система единиц неважна

> Fq:=(v2-v1)/dotprod(v2-v1,v2-v1)^(3/2);

> dotprod(v1,v2,'orthogonal');

> crossprod(v1,v2);

Уравнение 1: Ft перпендикулярен вектору v2

> eq1:=dotprod(Ft,v2,'orthogonal')=0;

Уравнение 2: Ft лежит в плоскости, образованной v1 и v2, т.е., он перпендикулярен векторному произведению v1 и v2.

> eq2:=dotprod(crossprod(v1,v2),Ft,'orthogonal')=0;

Уравнение 3: Ft - проекция Fq на направление Ft

> eq3:=dotprod(Fq,Ft,'orthogonal')=dotprod(Ft,Ft);

> so1:=solve({eq1,eq2,eq3},{ftx,fty,ftz});

> ftxf:=unapply(subs(so1[2],ftx),x1,y1,z1,x2,y2,z2):

> ftyf:=unapply(subs(so1[2],fty),x1,y1,z1,x2,y2,z2):

> ftzf:=unapply(subs(so1[2],ftz),x1,y1,z1,x2,y2,z2):

> df:=unapply(sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2),x1,y1,z1,x2,y2,z2):

Рассчитаем "траектории" зарядов, приняв, за Аристотелем, что скорость пропорциональна силе, приложенной к телу.

Заодно, рассчитаем и "потенциал" взаимодействия (также, в произвольных единицах, но энергия взаимодействия обратно

пропорциональна расстоянию, как и на самом деле)

> ct:=0.05:n:=5:m:=1500:

> x[1,0]:=0.:y[1,0]:=0.:z[1,0]:=1.:
x[2,0]:=cos(3.1415927/8.):y[2,0]:=sin(3.1415927/8.):z[2,0]:=0.:
x[3,0]:=0.:y[3,0]:=1.:z[3,0]:=0.:
x[4,0]:=-1.:y[4,0]:=0.:z[4,0]:=0.:
x[5,0]:=0.:y[5,0]:=-1.:z[5,0]:=0.:

> for t from 1 by 1 to m do
pot[t]:=0:
for i from 1 by 1 to n do
x[i,t] := x[i,t-1]:
y[i,t] := y[i,t-1]:
z[i,t] := z[i,t-1]:
fx[i,t]:=0:fy[i,t]:=0:fz[i,t]:=0:
for j from 1 by 1 to n do
if (i <> j) then
fx[i,t]:=fx[i,t]+ftxf(x[j,t-1],y[j,t-1],z[j,t-1],x[i,t-1],y[i,t-1],z[i,t-1]):
fy[i,t]:=fy[i,t]+ftyf(x[j,t-1],y[j,t-1],z[j,t-1],x[i,t-1],y[i,t-1],z[i,t-1]):
fz[i,t]:=fz[i,t]+ftzf(x[j,t-1],y[j,t-1],z[j,t-1],x[i,t-1],y[i,t-1],z[i,t-1]):
end if:
if (i < j) then
pot[t]:=pot[t]+1/df(x[j,t-1],y[j,t-1],z[j,t-1],x[i,t-1],y[i,t-1],z[i,t-1]):
end if:
end do:
x[i,t] := x[i,t] + ct * fx[i,t]:
y[i,t] := y[i,t] + ct * fy[i,t]:
z[i,t] := z[i,t] + ct * fz[i,t]:
end do:
end do:

> xf:=unapply(x[iv,floor(tv)],iv,tv):yf:=unapply(y[iv,floor(tv)],iv,tv):
zf:=unapply(z[iv,floor(tv)],iv,tv):

> potf:=unapply(pot[floor(tv)],tv):

> xff:=unapply(fx[iv,floor(tv)],iv,tv):yff:=unapply(fy[iv,floor(tv)],iv,tv):
zff:=unapply(fz[iv,floor(tv)],iv,tv):

Проекция траекторий на плоскость xz

> plot([seq([xf(ii,tt),zf(ii,tt),tt=1..m],ii=1..n)],style=point,symbol=circle,
color=[red,green,blue,black,yellow,violet]);

Проекция траекторий на плоскость xy

> plot([seq([xf(ii,tt),yf(ii,tt),tt=1..m],ii=1..n)],style=point,symbol=circle,
color=[red,green,blue,black,yellow,violet]);

Проекция Ft на плоскость xz

> plot([seq([xff(ii,tt),zff(ii,tt),tt=1..m],ii=1..n)],style=point,symbol=circle,
color=[red,green,blue,black,yellow,violet]);

Проекция Ft на плоскость xy

> plot([seq([xff(ii,tt),yff(ii,tt),tt=1..m],ii=1..n)],style=point,symbol=circle,
color=[red,green,blue,black,yellow,violet]);

Изменение потенциала с эволюцией системы

> plot([potf(tt)],tt=1..m);

Окончательное размещение зарядов

> me:=m:
ai := seq(arrow([x[ii,me],y[ii,me],z[ii,me]], shape=arrow),ii=1..n):
display(ai, scaling=CONSTRAINED, axes=NORMAL, labels=[x,y,z]);

>

> x[2,m]^2+y[2,m]^2+z[2,m]^2;