这里正在做什么？ （使用数学识别）

Question

我知道这本身与编程并不完全相关，但程序员是最重要的 也许所有人都会认识到这一点。

我有以下内容（X 和 Y 是数组，都有 3 个元素），但我无法识别（虽然它让我想起了一些事情，但不完全是！）这里正在做什么。它会给其他人敲响警钟吗？

我猜你可以忽略下面的部分；鞋面可能应该放弃它......但我仍然看不到它。

起初它让我想起 3d 空间中的线性插值......

  SUBROUTINE TRII(X,Y,XR,YR)
DIMENSION X(3),Y(3)

D=X(1)*(X(2)**2-X(3)**2)+
 >    X(2)*(X(3)**2-X(1)**2)+
 >    X(3)*(X(1)**2-X(2)**2)

D1=Y(1)*(X(2)*X(3)**2-X(3)*X(2)**2)+
 >     Y(2)*(X(3)*X(1)**2-X(1)*X(3)**2)+
 >     Y(3)*(X(1)*X(2)**2-X(2)*X(1)**2)

D2=Y(1)*(X(2)**2-X(3)**2)+
 >     Y(2)*(X(3)**2-X(1)**2)+
 >     Y(3)*(X(1)**2-X(2)**2)

D3=X(2)*(Y(3)-Y(1))+
 >     X(1)*(Y(2)-Y(3))+
 >     X(3)*(Y(1)-Y(2))

A=D1/D
B=D2/D
C=D3/D

YR=A+B*XR+C*XR**2

RETURN
END

  SUBROUTINE TRIM(X,Y,XR,YR,XM,YM)
DIMENSION X(3),Y(3)

D=X(1)*(X(2)**2-X(3)**2)+
 >    X(2)*(X(3)**2-X(1)**2)+
 >    X(3)*(X(1)**2-X(2)**2)

D1=Y(1)*(X(2)*X(3)**2-X(3)*X(2)**2)+
 >     Y(2)*(X(3)*X(1)**2-X(1)*X(3)**2)+
 >     Y(3)*(X(1)*X(2)**2-X(2)*X(1)**2)

D2=Y(1)*(X(2)**2-X(3)**2)+
 >     Y(2)*(X(3)**2-X(1)**2)+
 >     Y(3)*(X(1)**2-X(2)**2)

D3=X(2)*(Y(3)-Y(1))+
 >     X(1)*(Y(2)-Y(3))+
 >     X(3)*(Y(1)-Y(2))

A=D1/D
B=D2/D
C=D3/D

XR=-B/(2.*C)
YR=A+B*XR+C*XR**2

XM=XR
IF(XR.GT.X(1).OR.XR.LT.X(3))XM=X(1)
YM=A+B*XM+C*XM**2
IF(YM.LT.Y(1))XM=X(1)
IF(YM.LT.Y(1))YM=Y(1)

RETURN
END

“>”是一个延续符号。

Answer 1

这是一种求解线性方程组的方法，特别是cramers 规则。另请参阅萨鲁斯规则。之后，您似乎从中构造了一个二次方程。

Answer 2

我不确定这是什么语言，但很明显这是某种二次方程的求解器。 XR 和 YR 表达式是一个致命的泄露：

XR = -B / (2.*C)
YR = A + B*XR + C*XR**2

不知道 X(1..3) 和 Y(1.然而，.3) 表达式是，但是，不可能对 A/B/C 系数代表的内容推断太多。很多东西都使用二次方程——给定半径的圆面积、给定距离的光强度等等。需要更多上下文数据。

---

更新：OP表示，出于保密原因，他不能透露太多具体信息。不过，这里有一些提示：

• 子例程返回什么？稍后如何使用这些结果？这可能会带来更好的见解。

• 看起来 Y(1) 是该计算结果的某种神奇下界。请注意，如果 YM 小于 Y(1)，则 XM 和 YM 均设置为X(1) 和 Y(1)。

• “D”表达式看起来像这样，采用更自然的语法：

d = x1 * [x2^2 - x3^2] + x2 * [x3^2 - x1^2] + x3 * [x1^1 - x2^2]
  d1 = y1 * [x2*x3^2 - x3*x2^2] + y2 * [x3*x1^2 - x1*x3^2] + y3 * [x1*x2^2 - x1*x2^2]
  d2 = y1 * [x2^2 - x3^2] + y2 * [x3^2 - x1^2] + y3 * [x1^2 - x2^2]
  d3 = x2 * [y3 - y1] + x1 * [y2 - y3] * x3 * [y1 - y2]

• 这看起来非常像某种矩阵运算； D 几乎可以肯定是“行列式”。但还有其他事物具有相同的数学关系。

Answer 3

此代码表示一种对三个 2d 点的插值/二次曲线拟合，以及计算间隔本身内此类拟合二次曲线的最小值或最大值的方法。我猜想 TRII 代表三重（点）插值，TRIM 代表三重（点）最小值或最大值。

更精确地说，TRII 解决了这个问题：- 找到一条通过点 (x1,y1)、(x2,y2) 和 (x3,y3) 的二次曲线，形式为 Y=A+BX+C X^2 并计算点 XR 处二次方程的 Y 值并返回为 YR。这基本上是在三个 2d 点之间平滑插值的一种方法。它通常用于寻找一组离散数据点的最大值或最小值的更好近似值。

所有 D、D1、D2、D3 的内容都是为了求解矩阵方程：

(1 X1 X1^2) *(A) = (Y1)

(1 X2 X2^2) *(B) = (Y2)

(1 X3 X3^2) *(C) = (Y3)

使用其他评论之一中提到的 Cramers 规则，D 是矩阵行列式，D1、D2、D3 是辅因子。

TRIM 再次计算二次方程 Y=A+BX+CX^2，然后找到该二次方程的最大值/最小值 (XM, YM)。这是通过首先找到二次方程有转折点的点来完成的：如果 F(X)=A+BX+CX^2，F'(XR)=B+2*C *XR=0，或XR=-B/2*C，YR=A+BXR+CXR^2。然后有一些逻辑强制返回的 XM、YM 最小值或最大值位于特定范围内。

代码：

XM=XR
。
。
。
IF(YM.LT.Y(1))YM=Y(1)

有点奇怪，因为如果我们假设 GT 和 LT 分别表示大于和小于，那么我们需要假设 X3'<'X1，否则条件 (XR.GT.X(1).OR.XR.LT.X(3)) 很简单，XM、YM 设置为 X1、Y1。

因此 X3'<'X1 并且条件表示如果二次最大/最小值超出区间 (X1,X3)，则像以前一样将 (XM,YM) 设置为 (X1,Y1)。如果不是，则如果 Y1 高于 Y 中的最小/最大值，则再次将 (XM,YM) 设置为 (X1, Y1)。

很难理解这意味着什么，我怀疑代码可能是错误的！有什么想法吗？

伊万

Answer 4

D 是矩阵的行列式：

        | x(1) x(1)² 1 |
D = det | x(2) x(2)² 1 |
        | x(3) x(3)² 1 |

在 D1 中，最右边的列已替换为 Y：

         | x(1) x(1)² Y(1) |
D1 = det | x(2) x(2)² Y(2) |
         | x(3) x(3)² Y(3) |

在 D2 和 D3 中，它分别是第一列和第二列。现在是不是更容易认出来了？对我来说，看起来很像使用克莱默法则来求解线性方程。

编辑：更准确地说：（A，B，C）是系统的解：

A + x(1)*B + x(1)²*C = Y(1)
A + x(2)*B + x(2)²*C = Y(2)
A + x(3)*B + x(3)²*C = Y(3)

YR是二次方程解的平方（nb，不同的x！）：

C*x² + B*x + A = 0

我觉得现在这应该很明显了，但我不太明白......

Answer 5

代码运行如下

例程 TRII 将三个点 (x,y) 的坐标作为输入，并使用 拉格朗日插值。还采用坐标 XR 作为输入。在 YR 中返回插值抛物线 XR 处的值。
我猜例程的名称来自“TRI”（克罗地亚语“三”（点））和“I”插值。

例程 TRIM 也计算相同的抛物线，并返回区间{X(1),X(3)}中的函数。名称来自“TRI”和“M”（最小值）

（我“真的”执行了程序）>）

请注意，这是 FORTRAN 代码，参数是通过引用传递的，因此结果以相同的参数返回（非常奇怪！）

编辑

只是为了好玩，让我们运行 TRII

TRII[X_, Y_, XR_] := 
  Module[{D0, D1, D2, D3, A, B, C}, 
     D0 = X[[1]]*(X[[2]]^2 - X[[3]]^2) + 
          X[[2]]*(X[[3]]^2 - X[[1]]^2) + 
          X[[3]]*(X[[1]]^2 - X[[2]]^2);
     D1 = Y[[1]]*(X[[2]]*X[[3]]^2 - X[[3]]*X[[2]]^2) + 
          Y[[2]]*(X[[3]]*X[[1]]^2 - X[[1]]*X[[3]]^2) + 
          Y[[3]]*(X[[1]]*X[[2]]^2 - X[[2]]*X[[1]]^2);
     D2 = Y[[1]]*(X[[2]]^2 - X[[3]]^2) + 
          Y[[2]]*(X[[3]]^2 - X[[1]]^2) + 
          Y[[3]]*(X[[1]]^2 - X[[2]]^2);
     D3 = X[[2]]*(Y[[3]] - Y[[1]]) + 
          X[[1]]*(Y[[2]] - Y[[3]]) + 
          X[[3]]*(Y[[1]] - Y[[2]]);
   A = D1/D0;
   B = D2/D0;
   C = D3/D0;
   Return[A + B*XR + C*XR^2];];

X = RandomReal[1, 3];
Y = RandomReal[1, 3];
Show[Plot[TRII[X, Y, x], {x, 0, 1}], 
 ListPlot[Transpose[{X, Y}], PlotMarkers -> Automatic]]