前言
最近數值分析課上老師給出一道作業題，題目的内容為：
某地區為估計某礦物的儲量，在該地區内進行勘探，得到如下資料：
| 編号 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 |
| x坐标/km | ||||||||||
| y坐标/km | ||||||||||
| 礦物體 厚度H/m | 13.72 | 25.80 | 8.47 | 25.27 | 22.32 | 15.47 | 21.33 | 14.49 | 24.83 | 26.19 |
| 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | |
| 23.28 | 26.48 | 29.14 | 12.04 | 14.58 | 19.95 | 23.73 | 15.35 | 18.01 | 16.29 |
試估計出此地區内（$1<x<4,1<y<5$ ）該礦物的儲量。
老師還給了提示：礦物體的厚度$H$是坐标$x,y$的二進制函數，即面密度函數為：$H = H(x, y) $， 根據二重積分可知，所求礦物的儲量就是求二重積分$\iint\limits_{D}H\left ( x,y\right )dxdy$的值，其中$D=\left \{\left ( x,y\right )|1\leqslant x\leqslant 4,1\leqslant y\leqslant 5\right \}$。可考慮将二重積分化為累次積分，利用複合求積公式計算。（不相信我們的實力哈哈😉）大家看着這道題目會怎麼去做呢？
一、二進制拉格朗日插值
我首先想到的是用類似于一進制拉格朗日插值的二進制插值來對資料進行插值計算，得出一個形如$\sum_{i}^{m}\sum_{j}^{n}a_{ij}x^{i}y^{j}$的二進制函數，再對這個函數進行累次積分即可得出結果。根據一進制拉格朗日插值可推測二進制拉格朗日插值具有類似的形式：
其中，$f\left ( x_i,y_j \right )$為$\left ( x_i,y_j \right )$處表格中的值，$\sigma _i\left ( x \right )=\left ( x-x_0\right )\cdots \left ( x-x_{i-1}\right )\left ( x-x_{i+1}\right )\cdots \left ( x-x_m\right )$,$\omega _j\left ( y \right )=\left ( y-y_0\right )\cdots \left ( y-y_{j-1}\right )\left ( y-y_{j+1}\right )\cdots \left ( y-y_n\right )$。這裡我使用c#程式設計語言使用WinForm連接配接Mathematica來設計算法：
//二進制插值函數
public void BinaryInterpolation(string[] xdatastr, string[] ydatastr, string[] zdatastr)
{
string messagestr = "";
string pictpathstr = Directory.GetCurrentDirectory().ToString() + "\\images\\picture_" + imagenumber + ".png";
mathKernel1.CaptureMessages = true;
mathKernel1.CaptureGraphics = true;
mathKernel1.GraphicsHeight = pictureBox1.Height;
mathKernel1.GraphicsWidth = pictureBox1.Width;
mathKernel1.Compute("ExportString[SetPrecision[Expand[Simplify[{xdata = " + ArrayToTableFunction(xdatastr, 1) + ", ydata = " + ArrayToTableFunction(ydatastr, 1) + ", zdata = " + ArrayToTableFunction(zdatastr, ydatastr.Length) + ", Ln = 0, For[m = 1, m <= " + xdatastr.Length + ", m++, {\\[Sigma] = 1, \\[Sigma]k = 1, For[i = 1, i <= " + xdatastr.Length + ", i++, If[i != m, {\\[Sigma] = \\[Sigma] * (x - xdata[[i]]), \\[Sigma]k = \\[Sigma]k * (xdata[[m]] - xdata[[i]])}]], For[n = 1, n <= " + ydatastr.Length + ", n++, {\\[Omega] = 1, \\[Omega]k = 1, For[j = 1, j <= " + ydatastr.Length + ", j++, If[j != n, {\\[Omega] = \\[Omega] * (y - ydata[[j]]), \\[Omega]k = \\[Omega]k * (ydata[[n]] - ydata[[j]])}]], Ln = Ln + zdata[[m, n]]*\\[Sigma]/\\[Sigma]k*\\[Omega]/\\[Omega]k}]}], Ln}[[6]][[1]]]]," + significantdigits + "], \"MathML\"]");
result = mathKernel1.Result.ToString();
mathMLControl1.MC_loadXML(result);
foreach (string me in mathKernel1.Messages)
{
messagestr += me;
}
textBox4.Text = messagestr;
}
//将數組轉化為清單
public string ArrayToTableFunction(string[] data1, string[] data2)
{
StringBuilder outputstr = new StringBuilder();
outputstr.Append("{");
for (int i = 0; i < data1.Length; i++)
{
if (i == data1.Length - 1)
{
outputstr.Append("{" + data1[i] + "," + data2[i] + "}");
}
else
{
outputstr.Append("{" + data1[i] + "," + data2[i] + "},");
}
}
outputstr.Append("}");
return outputstr.ToString();
} 程式運作結果如下：
插值得出的函數為$-1.8415277x^3y^4 + 21.386389x^3y^3 - 83.882639x^3y^2 + 129.51277x^3y - 68.041667x^3 + 12.159166x^2y^4 - 140.705x^2y^3 + 548.72833x^2y^2 - 841.5375x^2y + 441.585x^2 - 21.029305xy^4 + 241.22527xy^3 - 927.47903xy^2 + 1397.153xy - 728.74333x + 5.8191667y^4 - 62.391667y^3 + 213.15083y^2 - 267.81833y + 146.47$
其圖像為：
對插值函數分别對$x$和$y$積分可得結果為252.193。
二、累次拉格朗日插值
做完二進制插值，我突然想到，既然積分有重積分和累次積分，那麼類似于累次積分，拉格朗日插值可否累次進行呢？帶着這個疑問，我開始進行了實驗：首先将表格寫成$x,y$坐标形式如下
| x/y | |||||
分别對每一行進行一進制插值，得出的插值函數對$y$進行積分，可以得到與$x$對應的積分結果，然後再對結果進行一進制插值并對$x$積分得出最終結果。利用程式實作：
//累次拉格朗日插值函數
public void LagrangePointSetsNIntegrationForm(string[] xdatastr, string[] ydatastr, string[] zdatastr)
{
string messagestr = "";
string[] firstintrgrate = new string[xdatastr.Length];
string[][] zdatameshgrid = new string[xdatastr.Length][];
mathKernel1.CaptureMessages = true;
for (int i = 0; i < xdatastr.Length; i++)
{
zdatameshgrid[i] = new string[ydatastr.Length];
}
for (int n = 0; n < zdatastr.Length; n++)
{
zdatameshgrid[n / ydatastr.Length][n % ydatastr.Length] = zdatastr[n];
}
for (int i = 0; i < xdatastr.Length; i++)
{
mathKernel1.Compute("Integrate[Expand[Simplify[{xdata =" + ListToArrayFunction(ydatastr) + ",ydata = " + ListToArrayFunction(zdatameshgrid[i]) + ", Ln = 0, For[i = 1, i <= " + ydatastr.Length + ",i++, {\\[Omega] = 1, \\[Omega]k = 1, For[j = 1, j <= " + ydatastr.Length + ", j++, If[j != i, \\[Omega] = \\[Omega] * (x - xdata[[j]])]], For[k = 1, k <= " + ydatastr.Length + ", k++, If[k != i, \\[Omega]k = \\[Omega]k * (xdata[[i]] -xdata[[k]])]], Ln = Ln + ydata[[i]]*\\[Omega]/\\[Omega]k}], Ln}[[5]]]],{x, " + ydatastr[0] + ", " + ydatastr[ydatastr.Length - 1] + "}]");
firstintrgrate[i] = mathKernel1.Result.ToString();
textBox6.AppendText(mathKernel1.Result.ToString() + " ");
}
for (int j = 0; j < xdatastr.Length; j++)
{
mathKernel1.Compute("Integrate[Expand[Simplify[{xdata =" + ListToArrayFunction(xdatastr) + ",ydata = " + ListToArrayFunction(firstintrgrate) + ", Ln = 0, For[i = 1, i <= " + xdatastr.Length + ",i++, {\\[Omega] = 1, \\[Omega]k = 1, For[j = 1, j <= " + xdatastr.Length + ", j++, If[j != i, \\[Omega] = \\[Omega] * (x - xdata[[j]])]], For[k = 1, k <= " + xdatastr.Length + ", k++, If[k != i, \\[Omega]k = \\[Omega]k * (xdata[[i]] -xdata[[k]])]], Ln = Ln + ydata[[i]]*\\[Omega]/\\[Omega]k}], Ln}[[5]]]],{x, " + xdatastr[0] + ", " + xdatastr[xdatastr.Length - 1] + "}]");
textBox5.Text = mathKernel1.Result.ToString();
}
foreach (string me in mathKernel1.Messages)
{
messagestr += me;
}
}
//将清單轉化為數組
public string ListToArrayFunction(string[] data)
{
StringBuilder outputstr = new StringBuilder();
outputstr.Append("{");
for (int i = 0; i < data.Length; i++)
{
if (i == data.Length - 1)
{
outputstr.Append(data[i]);
}
else
{
outputstr.Append(data[i] + ",");
}
}
outputstr.Append("}");
return outputstr.ToString();
} 結果顯示，累次進行的插值結果與二進制插值結果是相同的。那如果我改變插值順序，即先對$x$插值，再對$y$插值，結果會發生變化嗎？我将資料表進行轉置：
| y/x | ||||
輸入資料後運作程式：
從結果上看，雖然插值的順序不同，第一次的積分值不相同，但是第二次積分值确實相同的，為了避免偶然事件的影響，我已使用多組資料證明了這個結論的正确性。但遺憾的是，我暫時沒有找到理論依據，無法給出證明。
三、結論
拉格朗日插值和積分一樣，存在二維插值和累次插值，并且二者是相等的（可能需要滿足某些條件）。雖然目前缺少理論依據，但是對于數學模組化等的應用可以提供一個新方法。