關於二元插值問題的探討

前言

       最近數值分析課上老師給出一道作業題，題目的內容為：

    某地區為估計某礦物的儲量，在該地區內進行勘探，得到如下數據：

 表1  地區勘探數據表：

編號	01	02	03	04	05	06	07	08	09	10
x坐標/km	1	1	1	1	1	2	2	2	2	2
y坐標/km	1	2	3	4	5	1	2	3	4	5
礦物體 厚度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
x坐標/km	3	3	3	3	3	4	4	4	4	4
y坐標/km	1	2	3	4	5	1	2	3	4	5
礦物體 厚度H/m	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	1	2	3	4	5
1	13.72	25.80	8.47	25.27	22.32
2	15.47	21.33	14.49	24.83	26.19
3	23.28	26.48	29.14	12.04	14.58
4	19.95	23.73	15.35	18.01	16.29

分別對每一行進行一元插值，得出的插值函數對$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	1	2	3	4
1	13.72	15.47	23.28	19.95
2	25.80	21.33	26.48	23.73
3	8.47	14.49	29.14	15.35
4	25.27	24.83	12.04	18.01
5	22.32	26.19	14.58	16.29

輸入數據后運行程序：

       從結果上看，雖然插值的順序不同，第一次的積分值不相同，但是第二次積分值確實相同的，為了避免偶然事件的影響，我已使用多組數據證實了這個結論的正確性。但遺憾的是，我暫時沒有找到理論依據，無法給出證明。

三、結論

       拉格朗日插值和積分一樣，存在二維插值和累次插值，並且二者是相等的（可能需要滿足某些條件）。雖然目前缺少理論依據，但是對於數學建模等的應用可以提供一個新方法。