Kirchhoff积分公式验证
mathematica吧
全部回复
仅看楼主
level 6
冉qx 楼主
为了验证Kirchhoff积分公式,代码如下:
In[573]:= \[Lambda] = 10;
k = (2 \[Pi])/\[Lambda];
n = {x, y, z};
R = {1000, 0, 0};
\[Psi] = Exp[-I k Sqrt[(x - R[[1]])^2 + y^2 + z^2]]/Norm[R];
\[ScriptCapitalE] = ({1, -1, 0} Exp[(I k (x + y))/Sqrt[2]])/Sqrt[2];
t1 = 1/(4 Pi) (NSurfaceIntegrate[\[ScriptCapitalE] (
Grad[\[Psi],{x, y, z}] . n) - \[Psi] (\!\(
\*SubscriptBox[\(\[PartialD]\), \({{x, y, z}}\)]\[ScriptCapitalE]\) .
n), {x, y, z} \[Element] Sphere[],
Method -> "LocalAdaptive"] // Norm)
t2 = 1/(4 Pi) (NSurfaceIntegrate[
Grad[(\[Psi] n . \[ScriptCapitalE]),{x, y, z}] -
n (Grad[\[Psi],{x, y, z}] . \[ScriptCapitalE]), {x, y,
z} \[Element] Sphere[], Method -> "LocalAdaptive"] // Norm
t3 = -(1/(
4 Pi)) (NSurfaceIntegrate[((Div[\[ScriptCapitalE],{x, y,
z}]) n + \[ScriptCapitalE]\[Cross](Curl[
n,{x, y, z}])) \[Psi], {x, y, z} \[Element] Sphere[],
Method -> "LocalAdaptive"] // Norm)
Out[579]= 0.000129072
Out[580]= 0.000129072
Out[581]= 0.
按照图片的说法,t2,t3应该为零,但t2怎么不为零。图片如下:
2025年10月09日 06点10分 1
吧务
level 15
吃饱了撑的又靠AI辅助把这个问题研究了一遍,不得不承认现在AI是真好使,虽然到头来也没搞清楚你这是什么教材,但是对基尔霍夫积分定理的解释倒是很清楚。(但是,再强调一遍:
这绝对不是你提问时不写明背景信息的理由。)
首先贴上你那张截图里的各符号的含义(全是Gemini生成的,但我姑且检查了一遍,似乎没问题):
你用于验算公式的是个平面波:
带帽子的x和y(OverHat[x]和OverHat[y])是直角系的基向量,j是虚数单位。
然后是你代码的问题:
1. 你对形如 v·∇ 的记号的理解是错的,算子的乘法不可交换,v·∇ 和∇·v 不是一个东西。
2. 当被积函数是个矢量时,NSurfaceIntegrate 默认进行第二类曲面积分,而你这里需要的是对每一项进行单独积分。
2.1. 按2改过来之后会发现Method -> "LocalAdaptive"此时是不可用的,应该直接去掉。这个积分也不需要进行额外的选项调节。
3. 对于你这个算例,格林函数的分母部分准确地说应该是Sqrt[(x - R[[1]])^2 + y^2 + z^2],你似乎是利用了“R数值很大”这点进行近似处理(为了规避奇点?),这完全是多余的,NSurfaceIntegrate 能算好这个积分。
3.1. 但是,最好还是加个Method -> {Automatic, SymbolicProcessing -> 0}提升积分速度。
4. 如果你希望t1能够等于电场在观察点处的数值,观察点的位置就不能是球面外,因为对于平面波,闭曲面外任意点处用基尔霍夫积分公式算出来都是0。
4.1. 如果要算球内的情况,注意法向量的指向也要相应修改。
5. 截图中的“第二和第三项面积分为零”这句很有意思,也不知道是不是教材作者有意为之的:它也可以理解成“第二项和第三项之和为零”——如果是这样的话,那它就没问题,这个问题具体参考 黄晓伟, 盛新庆. 矢量基尔霍夫公式经典证明的漏洞与新的严格证明[J]. 物理学报, 2017, 66(16): 164201. DOI: 10.7498/aps.66.164201
总之,把上面这些内容全改过来之后:
2026年05月31日 04点05分 3
吧务
level 15
算了,把代码文本也贴了吧:
λ = 10;
k = (2 π)/λ;
n = -{x, y, z};
R = {1/2, 0, 0};
ψ = Exp[-I k Sqrt[(x - R[[1]])^2 + y^2 + z^2]]/
Sqrt[(x - R[[1]])^2 + y^2 + z^2];
ℰ = ({1, -1, 0} Exp[(I k (x + y))/Sqrt[2]])/Sqrt[2];
t1 = 1/(4 π)
NSurfaceIntegrate[#, {x, y, z} ∈ Sphere[],
Method -> {Automatic,
SymbolicProcessing ->
0}] & /@ (ℰ Grad[ψ,{x, y, z}] .
n - ψ D[ℰ, {{x, y, z}}] . n)
(* {0.689256 + 0.155678 I, -0.689256 - 0.155678 I, 0.} *)
ℰ /. {x -> R[[1]], y -> 0.}
(* {0.689731 + 0.155791 I, -0.689731 - 0.155791 I, 0} *)
t2 =
1/(4 π)
NSurfaceIntegrate[#, {x, y, z} ∈ Sphere[],
Method -> {Automatic,
SymbolicProcessing ->
0}] & /@ (Grad[(ψ n . ℰ),{x, y, z}] -
n Grad[ψ,{x, y, z}] . ℰ)
(* {-0.529117 + 0.3232 I, 0.529103 - 0.323199 I,
0.0000125048 + 1.7986*10^-6 I} *)
t3 = -(1/(4 π))
NSurfaceIntegrate[#, {x, y, z} ∈ Sphere[],
Method -> {Automatic, SymbolicProcessing -> 0}] & /@ ((Grad[
n,{x, y,
z}] . ℰ + ℰ\[Cross](Curl[
n,{x, y, z}])) ψ)
(* {0.52911 - 0.3232 I, -0.52911 + 0.3232 I, 0.} *)
t2 + t3
(* {-6.89144*10^-6 + 4.9005*10^-7 I, -6.89144*10^-6 +
4.9005*10^-7 I, 0.0000125048 + 1.7986*10^-6 I} *)
2026年05月31日 04点05分 4
1