格林公式应用实例

火力运用过程中的一个关键问题是计算多枚导弹打击复杂形状面目标有效毁伤面积。复杂形状面目标可以用任意多边形描述。有效毁伤面积等于多枚导弹联合覆盖区域与多边形目标区域重叠部分面积。此问题通常用数值算法计算。在有效毁伤区域边界上计算环路积分,可以用解析算法解决此问题,我们计算如下这样一个图形的面积。

计算步骤:

第一步:确定方向及端点

(1)确定方向:外逆内顺

(2)确定端点坐标:计算联合有效覆盖区域的每一段边界的端点坐标(交点坐标,圆心坐标,半径,可以通过电脑大致定坐标。)

假设各交点坐标分别为:

\(A(a_1,a_2), B(b_1,b_2), C(c_1,c_2), D(d_1,d_2), E(e_1,e_2), F(f_1,f_2), G(g_1,g_2), H(h_1,h_2), K(k_1,k_2)\)

圆点坐标分别为:

\(O_1(O_{11},O_{12}),O_2(O_{21},O_{22}), O_3(O_{31},O_{32})\)
假设圆的半径分别为:
\(R_1, R_2, R_3\)

第二步:计算两类曲线积分

用代码根据公式计算两类曲线积分的结果

沿线段曲线积分公式为

\(\frac{1}{2}\int_lxdy-ydx=\frac{x_iy_j-x_jy_i}2\)

x=c(1,1.1,3.5,2.5,3,3.7,4.2,3,2.8) #各交点横坐标
y=c(3.5,3,3.1,1.5,0.53,0.7,2.1,3.8,4) #各交点纵坐标
ox=c(1.9,3.6,3.2) #圆点横坐标
oy=c(3.5,3.2,1.52)#圆点纵坐标
#半径
r1=r2=r3=1
#
thetao1x=NULL
thetao1y=NULL
for (i in c(1,2,9)) {
  thetao1x=cbind(thetao1x,acos((x[i]-ox[1])/r1))
  thetao1y=cbind(thetao1y,acos((y[i]-oy[1])/r1))
}

thetao3x=NULL
thetao3y=NULL
for (j in c(4,5,6)) {
  thetao3x=cbind(thetao3x,acos((x[j]-ox[3])/r3))
  thetao3y=cbind(thetao3y,acos((y[j]-oy[3])/r3))
}
thetao2x=NULL
thetao2y=NULL
for (k in c(7,8)) {
  thetao2x=cbind(thetao2x,acos((x[i]-ox[2])/r2))
  thetao2y=cbind(thetao2y,acos((y[i]-oy[2])/r2))
}

#沿线段曲线积分
S11=(x[9]*y[1]-x[1]*y[9])/2
S12=(x[3]*y[2]-x[2]*y[3])/2
S13=(x[3]*y[4]-x[4]*y[3])/2

print(S11)
## [1] 2.9
print(S12)
## [1] 3.545
print(S13)
## [1] -1.25

沿弧段曲线积分公式为

\(\frac{1}{2}\int_lxdy-ydx=\frac{x_{O_i}(y_j-y_i)+y_{O_i}(x_i-x_j)+R_i^2(\theta_j-\theta_i)}2\)

#沿弧线曲线积分
S21=(ox[1]*(y[2]-y[1])+oy[1]*(x[2]-x[1])+r1^2*(thetao1x[2]-thetao1x[1]))/2
S22=(ox[3]*(y[4]-y[5])+oy[3]*(x[4]-x[5])+r3^2*(thetao3x[2]-thetao3x[1]))/2

S23=(ox[3]*(y[5]-y[6])+oy[3]*(x[5]-x[6])+r3^2*(thetao3x[3]-thetao3x[2]))/2
S24=(ox[2]*(y[7]-y[8])+oy[2]*(x[7]-x[8])+r3^2*(thetao2x[2]-thetao2x[1]))/2
S25=(ox[1]*(y[8]-y[9])+oy[1]*(x[8]-x[9])+r1^2*(thetao2x[2]-thetao1x[3]))/2
print(S21)
## [1] -0.3962371
print(S22)
## [1] 0.8849802
print(S23)
## [1] -1.166478
print(S24)
## [1] -1.14
print(S25)
## [1] 1.183532

第三步:累计求面积和

\(S_D=沿线段曲线积分+沿弧线曲线积分\)

#总面积
S=S11+S12+S13+S21+S22+S23+S24+S25
print(S)
## [1] 4.560797