Wolfram 语言

符号与数值微积分

求球面上的电荷分布

找出球面上使库仑势为最小的自由运动等电荷粒子的位置,这便是平衡时的电荷分布.

n 表示粒子数.

In[1]:=
Click for copyable input
n = 50;

{xi, yi, zi}^(th) 粒子的笛卡尔坐标.

In[2]:=
Click for copyable input
vars = Join[Array[x, n], Array[y, n], Array[z, n]];

目标是使库仑势最小.

In[3]:=
Click for copyable input
potential = Sum[((x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2)^-(1/ 2), {i, 1, n - 1}, {j, i + 1, n}];

由于粒子是在球面上,它们的坐标必须满足大小为 1 的约束.

In[4]:=
Click for copyable input
sphereconstr = Table[x[i]^2 + y[i]^2 + z[i]^2 == 1, {i, 1, n}];

使用球面坐标在球面上随机选择初始点.

In[5]:=
Click for copyable input
rpts = ConstantArray[1, n]; thetapts = RandomReal[{0, Pi}, n]; phipts = RandomReal[{-Pi, Pi}, n]; spherpts = Transpose[{rpts, thetapts, phipts}];

将初始点变换成笛卡尔坐标.

In[6]:=
Click for copyable input
cartpts = CoordinateTransform["Spherical" -> "Cartesian", spherpts];

重新排列初始点使其与变量顺序匹配.

In[7]:=
Click for copyable input
initpts = Flatten[Transpose[cartpts]];

以球面为约束将库仑势最小化.

In[8]:=
Click for copyable input
sol = FindMinimum[{potential, sphereconstr}, Thread[{vars, initpts}]];

从解当中提取出粒子的平衡位置.

In[9]:=
Click for copyable input
solpts = Table[{x[i], y[i], z[i]}, {i, 1, n}] /. sol[[2]];

图示结果.

In[10]:=
Click for copyable input
Show[ListPointPlot3D[solpts, PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}, {-1.1, 1.1}}, PlotStyle -> {{PointSize[.03], Blue}}, AspectRatio -> 1, BoxRatios -> 1, PlotLabel -> "Particle Distribution"], Graphics3D[{Opacity[.5], Sphere[]}]]
Out[10]=

相关范例

de en es fr ja ko pt-br ru