[Numerical_analysis]the_third_project

一、题目

关于x，y，t，u，v，w的下列方程组

$$0.5cost+u+v+w-x=2.67$$
$$t+0.5sinu+v+w-y=1.07$$
$$0.5t+u+cosv+w-x=3.74$$
$$t+0.5u+v+sinw-y=0.79$$
以及关于z，t，u的二维数表确定了一个二元函数z=f(x，y)。
1．试用数值方法求出f(x，y)在区域D={(x,y)︱0≤x≤0.8,0.5≤y≤1.5}上的一个近似表达式 $$p(x,y)=\sum_{s=0}^kx^ry^s$$ 要求p(x,y)一最小的k值达到以下的精度
$$\sigma=\sum_{i=0}^{10}\sum_{j=0}^{20}[f(x_i-y_j)-p(x_i,y_j)]^2\leq10^{-7}$$ 其中xi=0.08i，yj=0.5+0.05j。
2．计算f(xi*,yj*),p(xi*,yj*)(i=1,2,…,8;j=1,2,…,5)的值，以观察p(x,y)逼近f(x,y)的效果，其中xi*=0.1i,yj*=0.5+0.2j。

二、算法设计方案

1.将$x_i^*=0.08i(i=0,1,2,...10;)$
和$y_j^*=0.5+0.05j(j=0,1,2,...20)$代入非线性方程组中，用牛顿法解出$t_i$和$u_j$；
2.以采取分片二次插值，选择（m，n）满足

$$t_m-h/2<t_i\leq t_m+h/2,2\leq m\leq 3$$ $$u_n-\tau/2<u_j\leq u_n+\tau/2,2\leq n\leq3$$ 如果 $t_i\leq t_1+h/2$或 $t_i> t_4-h/2$，则m=1或4；如果 $u_j\leq u_1+\tau/2$或 $u_j>u_4-\tau/2$，则n=1或n=4。选择 $(t_k,u_r)(k=m-1,m,m+1;r=n-1,n,n+1)$为插值节点，相应的Lagrange形式的插值多项式为 $$p_{22}(t,u)=\sum_{k=m-1}^{m+1}\sum_{r=n-1}^{n+1}l_k(t)\hat l_rf(t_k,u_r)$$ 其中 $$l_k(t)=\prod_{w=m-1 ; w\neq r}^{m+1}{t-t_w\over t_k-t_w}(k=m-1, m, m+1)$$ $$l_r(u)=\prod_{w=n-1;w\neq r}^{n+1}{y-y_w \over y_r-y_w} (r=n-1, n, n+1)$$
并将 $t_I$和 $u_j$代入 $p_{22}(t,u)$，便得到了数表 $x_i,y_j,f(x_i,y_j)$。
3.进行曲面拟和系数矩阵 $\mathbf C=[c_rs^*]$， $\mathbf C=(\mathbf B^\mathbf T\mathbf B)^{-1}\mathbf B^\mathbf T\mathbf U\mathbf G(\mathbf G^\mathbf T\mathbf G)^{-1}$
其中 $$\mathbf B=[\phi_c(x_i)]= \left[ \begin{matrix} 1 &x_0&L&x_k^k\\ 1&x_1&L&x_1^k\\ M&M&M&M\\ 1&x_{10}&L&x_{10}^k \end{matrix} \right]$$
，
$$\mathbf G=[\psi _s(y_j)]= \left[ \begin{matrix} 1 &y_0&L&y_k^k\\ 1&y_1&L&y_1^k\\ M&M&M&M\\ 1&y_{10}&L&y_{10}^k \end{matrix} \right]$$
$$\mathbf U=[f(x_i,y_j)]$$

k从0逐渐增大，直到$\sigma \leq 10^{-7}$，便得到了要求精度的系数$c_rs$。

三、运算结果

得到插值节点：
x(0)=0,y(0)=0,z= 4.46504018481e-01;x(0)=0,y(1)=1,z= 3.24683262928e-01;
x(0)=0,y(2)=1,z= 2.10159686683e-01;x(0)=0,y(3)=1,z= 1.03043608316e-01;
x(0)=0,y(4)=1,z= 3.40189556268e-03;x(0)=0,y(5)=1,z=-8.87358136380e-02;
x(0)=0,y(6)=1,z=-1.73371632750e-01;x(0)=0,y(7)=1,z=-2.50534611467e-01;
x(0)=0,y(8)=1,z=-3.20276506388e-01;x(0)=0,y(9)=1,z=-3.82668066110e-01;
x(0)=0,y(10)=1,z=-4.37795766738e-01;x(0)=0,y(11)=1,z=-4.85758941444e-01;
x(0)=0,y(12)=1,z=-5.26667254884e-01;x(0)=0,y(13)=1,z=-5.60638479797e-01;
x(0)=0,y(14)=1,z=-5.87796538768e-01;x(0)=0,y(15)=1,z=-6.08269779090e-01;
x(0)=0,y(16)=1,z=-6.22189452876e-01;x(0)=0,y(17)=1,z=-6.29688378186e-01;
x(0)=0,y(18)=1,z=-6.30899760003e-01;x(0)=0,y(19)=1,z=-6.25956152545e-01;
x(0)=0,y(20)=2,z=-6.14988546609e-01;
x(1)=0,y(0)=0,z= 6.38015226511e-01;x(1)=0,y(1)=1,z= 5.06611755147e-01;
x(1)=0,y(2)=1,z= 3.82176369277e-01;x(1)=0,y(3)=1,z= 2.64863491154e-01;
x(1)=0,y(4)=1,z= 1.54780200285e-01;x(1)=0,y(5)=1,z= 5.19926834909e-02;
x(1)=0,y(6)=1,z=-4.34680402049e-02;x(1)=0,y(7)=1,z=-1.31601056789e-01;
x(1)=0,y(8)=1,z=-2.12431088309e-01;x(1)=0,y(9)=1,z=-2.86004551058e-01;
x(1)=0,y(10)=1,z=-3.52386078979e-01;x(1)=0,y(11)=1,z=-4.11655456522e-01;
x(1)=0,y(12)=1,z=-4.63904911519e-01;x(1)=0,y(13)=1,z=-5.09236724701e-01;
x(1)=0,y(14)=1,z=-5.47761117962e-01;x(1)=0,y(15)=1,z=-5.79594388339e-01;
x(1)=0,y(16)=1,z=-6.04857258890e-01;x(1)=0,y(17)=1,z=-6.23673421332e-01;
x(1)=0,y(18)=1,z=-6.36168248413e-01;x(1)=0,y(19)=1,z=-6.42467656690e-01;
x(1)=0,y(20)=2,z=-6.42697102700e-01;
x(2)=0,y(0)=0,z= 8.40081395767e-01;x(2)=0,y(1)=1,z= 6.99764165673e-01;
x(2)=0,y(2)=1,z= 5.66061442352e-01;x(2)=0,y(3)=1,z= 4.39171608118e-01;
x(2)=0,y(4)=1,z= 3.19242138041e-01;x(2)=0,y(5)=1,z= 2.06376192387e-01;
x(2)=0,y(6)=1,z= 1.00638523891e-01;x(2)=0,y(7)=1,z= 2.06074006784e-03;
x(2)=0,y(8)=1,z=-8.93540247670e-02;x(2)=0,y(9)=1,z=-1.73626968865e-01;
x(2)=0,y(10)=1,z=-2.50799956160e-01;x(2)=0,y(11)=1,z=-3.20932269445e-01;
x(2)=0,y(12)=1,z=-3.84097735005e-01;x(2)=0,y(13)=1,z=-4.40382175418e-01;
x(2)=0,y(14)=1,z=-4.89881152313e-01;x(2)=0,y(15)=1,z=-5.32697965534e-01;
x(2)=0,y(16)=1,z=-5.68941879292e-01;x(2)=0,y(17)=1,z=-5.98726549515e-01;
x(2)=0,y(18)=1,z=-6.22168629750e-01;x(2)=0,y(19)=1,z=-6.39386535697e-01;
x(2)=0,y(20)=2,z=-6.50499350788e-01;
x(3)=0,y(0)=0,z= 1.05151509180e+00;x(3)=0,y(1)=1,z= 9.02927430831e-01;
x(3)=0,y(2)=1,z= 7.60580266860e-01;x(3)=0,y(3)=1,z= 6.24715198146e-01;
x(3)=0,y(4)=1,z= 4.95519756001e-01;x(3)=0,y(5)=1,z= 3.73134042775e-01;
x(3)=0,y(6)=1,z= 2.57656748872e-01;x(3)=0,y(7)=1,z= 1.49150559410e-01;
x(3)=0,y(8)=1,z= 4.76469867734e-02;x(3)=0,y(9)=1,z=-4.68493232015e-02;
x(3)=0,y(10)=1,z=-1.34356760385e-01;x(3)=0,y(11)=1,z=-2.14913344927e-01;
x(3)=0,y(12)=1,z=-2.88573700635e-01;x(3)=0,y(13)=1,z=-3.55406364786e-01;
x(3)=0,y(14)=1,z=-4.15491396489e-01;x(3)=0,y(15)=1,z=-4.68918249969e-01;
x(3)=0,y(16)=1,z=-5.15783883125e-01;x(3)=0,y(17)=1,z=-5.56191075200e-01;
x(3)=0,y(18)=1,z=-5.90246930563e-01;x(3)=0,y(19)=1,z=-6.18061548241e-01;
x(3)=0,y(20)=2,z=-6.39746839258e-01;
x(4)=0,y(0)=0,z= 1.27124675148e+00;x(4)=0,y(1)=1,z= 1.11500201815e+00;
x(4)=0,y(2)=1,z= 9.64607727216e-01;x(4)=0,y(3)=1,z= 8.20347369475e-01;
x(4)=0,y(4)=1,z= 6.82447678179e-01;x(4)=0,y(5)=1,z= 5.51085208598e-01;
x(4)=0,y(6)=1,z= 4.26392385902e-01;x(4)=0,y(7)=1,z= 3.08462995633e-01;
x(4)=0,y(8)=1,z= 1.97357129692e-01;x(4)=0,y(9)=1,z= 9.31056208594e-02;
x(4)=0,y(10)=1,z=-4.28599223403e-03;x(4)=0,y(11)=1,z=-9.48339252969e-02;
x(4)=0,y(12)=1,z=-1.78572990364e-01;x(4)=0,y(13)=1,z=-2.55553779055e-01;
x(4)=0,y(14)=1,z=-3.25840150158e-01;x(4)=0,y(15)=1,z=-3.89506988363e-01;
x(4)=0,y(16)=1,z=-4.46638204599e-01;x(4)=0,y(17)=1,z=-4.97324951768e-01;
x(4)=0,y(18)=1,z=-5.41664032699e-01;x(4)=0,y(19)=1,z=-5.79756479795e-01;
x(4)=0,y(20)=2,z=-6.11706288148e-01;
x(5)=0,y(0)=0,z= 1.49832105248e+00;x(5)=0,y(1)=1,z= 1.33499863207e+00;
x(5)=0,y(2)=1,z= 1.17712512374e+00;x(5)=0,y(3)=1,z= 1.02502405502e+00;
x(5)=0,y(4)=1,z= 8.78960023174e-01;x(5)=0,y(5)=1,z= 7.39145108704e-01;
x(5)=0,y(6)=1,z= 6.05744871487e-01;x(5)=0,y(7)=1,z= 4.78883861067e-01;
x(5)=0,y(8)=1,z= 3.58650625882e-01;x(5)=0,y(9)=1,z= 2.45102236196e-01;
x(5)=0,y(10)=1,z= 1.38268350928e-01;x(5)=0,y(11)=1,z= 3.81548654070e-02;
x(5)=0,y(12)=1,z=-5.52528211681e-02;x(5)=0,y(13)=1,z=-1.41986880814e-01;
x(5)=0,y(14)=1,z=-2.22094439096e-01;x(5)=0,y(15)=1,z=-2.95635232460e-01;
x(5)=0,y(16)=1,z=-3.62679511503e-01;x(5)=0,y(17)=1,z=-4.23306164224e-01;
x(5)=0,y(18)=1,z=-4.77601036132e-01;x(5)=0,y(19)=1,z=-5.25655426667e-01;
x(5)=0,y(20)=2,z=-5.67564743655e-01;
x(6)=0,y(0)=0,z= 1.73189274038e+00;x(6)=0,y(1)=1,z= 1.56203457721e+00;
x(6)=0,y(2)=1,z= 1.39721691821e+00;x(6)=0,y(3)=1,z= 1.23780100674e+00;
x(6)=0,y(4)=1,z= 1.08408753268e+00;x(6)=0,y(5)=1,z= 9.36322772315e-01;
x(6)=0,y(6)=1,z= 7.94704449054e-01;x(6)=0,y(7)=1,z= 6.59387198028e-01;
x(6)=0,y(8)=1,z= 5.30487586840e-01;x(6)=0,y(9)=1,z= 4.08088685454e-01;
x(6)=0,y(10)=1,z= 2.92244201230e-01;x(6)=0,y(11)=1,z= 1.82982206854e-01;
x(6)=0,y(12)=1,z= 8.03084940354e-02;x(6)=0,y(13)=1,z=-1.57904130516e-02;
x(6)=0,y(14)=1,z=-1.05344551621e-01;x(6)=0,y(15)=1,z=-1.88398090610e-01;
x(6)=0,y(16)=1,z=-2.65007149319e-01;x(6)=0,y(17)=1,z=-3.35237838904e-01;
x(6)=0,y(18)=1,z=-3.99164503887e-01;x(6)=0,y(19)=1,z=-4.56868143302e-01;
x(6)=0,y(20)=2,z=-5.08434993278e-01;
x(7)=1,y(0)=0,z= 1.97122178640e+00;x(7)=1,y(1)=1,z= 1.79532959950e+00;
x(7)=1,y(2)=1,z= 1.62406711323e+00;x(7)=1,y(3)=1,z= 1.45783058271e+00;
x(7)=1,y(4)=1,z= 1.29695464975e+00;x(7)=1,y(5)=1,z= 1.14171810545e+00;
x(7)=1,y(6)=1,z= 9.92349533324e-01;x(7)=1,y(7)=1,z= 8.49032663329e-01;
x(7)=1,y(8)=1,z= 7.11911352264e-01;x(7)=1,y(9)=1,z= 5.81094158922e-01;
x(7)=1,y(10)=1,z= 4.56658513233e-01;x(7)=1,y(11)=1,z= 3.38654496139e-01;
x(7)=1,y(12)=1,z= 2.27108255770e-01;x(7)=1,y(13)=1,z= 1.22025089193e-01;
x(7)=1,y(14)=1,z= 2.33922196376e-02;x(7)=1,y(15)=1,z=-6.88187019710e-02;
x(7)=1,y(16)=1,z=-1.54649344213e-01;x(7)=1,y(17)=1,z=-2.34152666459e-01;
x(7)=1,y(18)=1,z=-3.07391091913e-01;x(7)=1,y(19)=1,z=-3.74434862348e-01;
x(7)=1,y(20)=2,z=-4.35360556536e-01;
x(8)=1,y(0)=0,z= 2.21566786369e+00;x(8)=1,y(1)=1,z= 2.03420113361e+00;
x(8)=1,y(2)=1,z= 1.85695514362e+00;x(8)=1,y(3)=1,z= 1.68435816416e+00;
x(8)=1,y(4)=1,z= 1.51677635240e+00;x(8)=1,y(5)=1,z= 1.35451904115e+00;
x(8)=1,y(6)=1,z= 1.19784408667e+00;x(8)=1,y(7)=1,z= 1.04696304942e+00;
x(8)=1,y(8)=1,z= 9.02046083802e-01;x(8)=1,y(9)=1,z= 7.63226477663e-01;
x(8)=1,y(10)=1,z= 6.30604821954e-01;x(8)=1,y(11)=1,z= 5.04252814597e-01;
x(8)=1,y(12)=1,z= 3.84216715546e-01;x(8)=1,y(13)=1,z= 2.70520476641e-01;
x(8)=1,y(14)=1,z= 1.63168572400e-01;x(8)=1,y(15)=1,z= 6.21485581168e-02;
x(8)=1,y(16)=1,z=-3.25666193968e-02;x(8)=1,y(17)=1,z=-1.21016534844e-01;
x(8)=1,y(18)=1,z=-2.03251399623e-01;x(8)=1,y(19)=1,z=-2.79330359558e-01;
x(8)=1,y(20)=2,z=-3.49319957540e-01;
x(9)=1,y(0)=0,z= 2.46468422266e+00;x(9)=1,y(1)=1,z= 2.27805897940e+00;
x(9)=1,y(2)=1,z= 2.09525125084e+00;x(9)=1,y(3)=1,z= 1.91671812800e+00;
x(9)=1,y(4)=1,z= 1.74285462878e+00;x(9)=1,y(5)=1,z= 1.57399842733e+00;
x(9)=1,y(6)=1,z= 1.41043483523e+00;x(9)=1,y(7)=1,z= 1.25240175061e+00;
x(9)=1,y(8)=1,z= 1.10009440963e+00;x(9)=1,y(9)=1,z= 9.53669851261e-01;
x(9)=1,y(10)=1,z= 8.13251055249e-01;x(9)=1,y(11)=1,z= 6.78930742966e-01;
x(9)=1,y(12)=1,z= 5.50774848504e-01;x(9)=1,y(13)=1,z= 4.28825676973e-01;
x(9)=1,y(14)=1,z= 3.13104771740e-01;x(9)=1,y(15)=1,z= 2.03615514033e-01;
x(9)=1,y(16)=1,z= 1.00345478241e-01;x(9)=1,y(17)=1,z= 3.26856518657e-03;
x(9)=1,y(18)=1,z=-8.76530659133e-02;x(9)=1,y(19)=1,z=-1.72467247819e-01;
x(9)=1,y(20)=2,z=-2.51230220752e-01;
x(10)=1,y(0)=0,z= 2.71781110947e+00;x(10)=1,y(1)=1,z= 2.52639950126e+00;
x(10)=1,y(2)=1,z= 2.33841138686e+00;x(10)=1,y(3)=1,z= 2.15432937728e+00;
x(10)=1,y(4)=1,z= 1.97457455665e+00;x(10)=1,y(5)=1,z= 1.79951057910e+00;
x(10)=1,y(6)=1,z= 1.62944822055e+00;x(10)=1,y(7)=1,z= 1.46465004375e+00;
x(10)=1,y(8)=1,z= 1.30533496765e+00;x(10)=1,y(9)=1,z= 1.15168262131e+00;
x(10)=1,y(10)=1,z= 1.00383741991e+00;x(10)=1,y(11)=1,z= 8.61912337228e-01;
x(10)=1,y(12)=1,z= 7.25992371111e-01;x(10)=1,y(13)=1,z= 5.96137711520e-01;
x(10)=1,y(14)=1,z= 4.72386627914e-01;x(10)=1,y(15)=1,z= 3.54758095898e-01;
x(10)=1,y(16)=1,z= 2.43254184181e-01;x(10)=1,y(17)=1,z= 1.37862222525e-01;
x(10)=1,y(18)=1,z= 3.85567703264e-02;x(10)=1,y(19)=1,z=-5.46985959345e-02;
x(10)=1,y(20)=2,z=-1.41949659709e-01;

k=0; Sigma=1.44288077184e+02;
k=1; Sigma=3.22090897364e+00;
k=2; Sigma=4.65996003327e-03;
k=3; Sigma=1.72117537914e-04;
k=4; Sigma=3.30953429925e-06;
k=5; Sigma=2.54137804282e-08;
满足要求的最小的k=5; 此时误差sigma=2.54137804282e-08;

拟合得到的多项式系数如下：
c(0,0)= 2.02123051866e+00; c(0,1)=-3.66842679642e+00; c(0,2)= 7.09248619309e-01; 
c(0,3)= 8.48605406376e-01; c(0,4)=-4.15897435997e-01; c(0,5)= 6.74319952231e-02; 
c(1,0)= 3.19190900891e+00; c(1,1)=-7.41110377179e-01; c(1,2)=-2.69712460103e+00; 
c(1,3)= 1.63118340107e+00; c(1,4)=-4.84719972985e-01; c(1,5)= 6.06142812009e-02; 
c(2,0)= 2.56889811764e-01; c(2,1)= 1.57991867306e+00; c(2,2)=-4.63408148115e-01; 
c(2,3)=-8.13443916114e-02; c(2,4)= 1.02094237837e-01; c(2,5)=-2.10152316721e-02; 
c(3,0)=-2.69260334105e-01; c(3,1)=-7.30247662788e-01; c(3,2)= 1.07614508077e+00; 
c(3,3)=-8.07012812330e-01; c(3,4)= 3.02872801715e-01; c(3,5)=-4.59726348959e-02; 
c(4,0)= 2.17459751388e-01; c(4,1)=-1.78372377248e-01; c(4,2)=-7.24058018635e-02; 
c(4,3)= 2.43330466155e-01; c(4,4)=-1.41334735658e-01; c(4,5)= 2.65102404470e-02; 
c(5,0)=-5.59032636069e-02; c(5,1)= 1.43199240644e-01; c(5,2)=-1.36270369467e-01; 
c(5,3)= 4.07196217619e-02; c(5,4)= 3.77503209727e-03; c(5,5)=-2.66770111942e-03; 

计算拟合值和实际值：
x*(1)=0,y*(1)=1,f(x*,y*)=  1.94720407918e-01,p(x*,y*)=  1.94730358437e-01;
x*(1)=0,y*(2)=1,f(x*,y*)= -1.83037079189e-01,p(x*,y*)= -1.83041836594e-01;
x*(1)=0,y*(3)=1,f(x*,y*)= -4.45497646915e-01,p(x*,y*)= -4.45500037580e-01;
x*(1)=0,y*(4)=1,f(x*,y*)= -5.97566707641e-01,p(x*,y*)= -5.97558849594e-01;
x*(1)=0,y*(5)=2,f(x*,y*)= -6.46459593901e-01,p(x*,y*)= -6.46446099038e-01;
x*(2)=0,y*(1)=1,f(x*,y*)=  4.05979189288e-01,p(x*,y*)=  4.05989540581e-01;
x*(2)=0,y*(2)=1,f(x*,y*)= -2.25159583746e-02,p(x*,y*)= -2.25211142351e-02;
x*(2)=0,y*(3)=1,f(x*,y*)= -3.38220816040e-01,p(x*,y*)= -3.38224018838e-01;
x*(2)=0,y*(4)=1,f(x*,y*)= -5.44437831522e-01,p(x*,y*)= -5.44430444430e-01;
x*(2)=0,y*(5)=2,f(x*,y*)= -6.47361338568e-01,p(x*,y*)= -6.47348000459e-01;
x*(3)=0,y*(1)=1,f(x*,y*)=  6.34777195151e-01,p(x*,y*)=  6.34787453510e-01;
x*(3)=0,y*(2)=1,f(x*,y*)=  1.58801168839e-01,p(x*,y*)=  1.58796297495e-01;
x*(3)=0,y*(3)=1,f(x*,y*)= -2.07365694171e-01,p(x*,y*)= -2.07368577378e-01;
x*(3)=0,y*(4)=1,f(x*,y*)= -4.65357906898e-01,p(x*,y*)= -4.65349918045e-01;
x*(3)=0,y*(5)=2,f(x*,y*)= -6.20270953075e-01,p(x*,y*)= -6.20257130500e-01;
x*(4)=0,y*(1)=1,f(x*,y*)=  8.78960023174e-01,p(x*,y*)=  8.78969865697e-01;
x*(4)=0,y*(2)=1,f(x*,y*)=  3.58650625882e-01,p(x*,y*)=  3.58646044354e-01;
x*(4)=0,y*(3)=1,f(x*,y*)= -5.52528211681e-02,p(x*,y*)= -5.52554347839e-02;
x*(4)=0,y*(4)=1,f(x*,y*)= -3.62679511503e-01,p(x*,y*)= -3.62671059281e-01;
x*(4)=0,y*(5)=2,f(x*,y*)= -5.67564743655e-01,p(x*,y*)= -5.67550576839e-01;
x*(5)=0,y*(1)=1,f(x*,y*)=  1.13661091016e+00,p(x*,y*)=  1.13662035338e+00;
x*(5)=0,y*(2)=1,f(x*,y*)=  5.74980340948e-01,p(x*,y*)=  5.74975843163e-01;
x*(5)=0,y*(3)=1,f(x*,y*)=  1.15992376792e-01,p(x*,y*)=  1.15989322501e-01;
x*(5)=0,y*(4)=1,f(x*,y*)= -2.38568304012e-01,p(x*,y*)= -2.38560417403e-01;
x*(5)=0,y*(5)=2,f(x*,y*)= -4.91434393656e-01,p(x*,y*)= -4.91420897528e-01;
x*(6)=1,y*(1)=1,f(x*,y*)=  1.40604179891e+00,p(x*,y*)=  1.40605068686e+00;
x*(6)=1,y*(2)=1,f(x*,y*)=  8.05941494063e-01,p(x*,y*)=  8.05937301650e-01;
x*(6)=1,y*(3)=1,f(x*,y*)=  3.04429221045e-01,p(x*,y*)=  3.04425831855e-01;
x*(6)=1,y*(4)=1,f(x*,y*)= -9.50161300996e-02,p(x*,y*)= -9.50089460349e-02;
x*(6)=1,y*(5)=2,f(x*,y*)= -3.93902307746e-01,p(x*,y*)= -3.93889837910e-01;
x*(7)=1,y*(1)=1,f(x*,y*)=  1.68578351531e+00,p(x*,y*)=  1.68579121679e+00;
x*(7)=1,y*(2)=1,f(x*,y*)=  1.04988115306e+00,p(x*,y*)=  1.04987773830e+00;
x*(7)=1,y*(3)=1,f(x*,y*)=  5.08293783940e-01,p(x*,y*)=  5.08291043158e-01;
x*(7)=1,y*(4)=1,f(x*,y*)=  6.61487967065e-02,p(x*,y*)=  6.61563528951e-02;
x*(7)=1,y*(5)=2,f(x*,y*)= -2.76834341778e-01,p(x*,y*)= -2.76822047537e-01;
x*(8)=1,y*(1)=1,f(x*,y*)=  1.97457455665e+00,p(x*,y*)=  1.97458126046e+00;
x*(8)=1,y*(2)=1,f(x*,y*)=  1.30533496765e+00,p(x*,y*)=  1.30533200220e+00;
x*(8)=1,y*(3)=1,f(x*,y*)=  7.25992371111e-01,p(x*,y*)=  7.25989307481e-01;
x*(8)=1,y*(4)=1,f(x*,y*)=  2.43254184181e-01,p(x*,y*)=  2.43260785216e-01;
x*(8)=1,y*(5)=2,f(x*,y*)= -1.41949659709e-01,p(x*,y*)= -1.41938797091e-01;

四、源程序

#include<math.h>
#include<iostream>
#include<stdio.h>
using namespaced std;
void DoLittle(double a[n][n],double b[n])
    {
    int i,j,k,t,i_k,m[n];
    double u[n][n],l[n][n],s[n],y[n];
    double u_x,l_u,max,temp;
    for(k=1;k<=n;k++)
    { 
        for(i=k;i<=n;i++)
        {   
            l_u=0;
            for(t=1;t<=k-1;t++)  
                l_u=l_u+l[i-1][t-1]*u[t-1][k-1];
            s[i-1]=a[i-1][k-1]-l_u;
        }
        max=fabs(s[k-1]);
        i_k=k;
        m[k-1]=k;
        for(i=k;i<=n;i++)
        {
            if(fabs(s[i-1])>max)
            {
                max=fabs(s[i-1]);
                i_k=i;
                m[k-1]=i_k;
            }
        }
        if(i_k!=k)
        {
            for(t=1;t<=k-1;t++)
            {   
                temp=l[k-1][t-1];
                l[k-1][t-1]=l[i_k-1][t-1];
                l[i_k-1][t-1]=temp;
            }
            for(t=k;t<=n;t++)
            {           
                temp=a[k-1][t-1];
                a[k-1][t-1]=a[i_k-1][t-1];
                a[i_k-1][t-1]=temp;
            }
            temp=s[k-1];
            s[k-1]=s[i_k-1];
            s[i_k-1]=temp;
        }
        u[k-1][k-1]=s[k-1];
        for(j=k+1;j<=n;j++)
        {       
            l_u=0;
            for(t=1;t<=k-1;t++) 
                l_u=l_u+l[k-1][t-1]*u[t-1][j-1];
            u[k-1][j-1]=a[k-1][j-1]-l_u;
        }
        for(i=k+1;i<=n;i++)  
            l[i-1][k-1]=s[i-1]/u[k-1][k-1];     
    }
    for(k=1;k<=n-1;k++) 
    {
        temp=b[k-1];
        b[k-1]=b[m[k-1]-1];
        b[m[k-1]-1]=temp;
    }
    y[0]=b[0];
    for(i=2;i<=n;i++)
    {
        l_u=0;
        for(t=1;t<=i-1;t++) 
            l_u=l_u+l[i-1][t-1]*y[t-1]; 
        y[i-1]=b[i-1]-l_u;
    }
    x_[n-1]=y[n-1]/u[n-1][n-1];
    for(i=n-1;i>=1;i--)
    {
        u_x=0;
        for(t=i+1;t<=n;t++)
            u_x=u_x+u[i-1][t-1]*x_[t-1];
        x_[i-1]=(y[i-1]-u_x)/u[i-1][i-1];
    }
}

double MAX(double a[n])//求数组中的最大值
{
    int i;
    double max;
    max=fabs(a[0]);
    for(i=0;i<n;i++)
        if(fabs(a[i])>max)
            max=fabs(a[i]);     
    return(max);
}

void NEWTON(double x_x,double y_y)
{ 
        double F[n],F_1[n][n];
    int i,j;
    for(i=0;i<n;i++)
    {
        x_[i]=x[i]=1;
    }
    for(i=0;i<2000;i++)
    {
        F[0]=(2.67-0.5*cos(x[0])-x[1]-x[2]-x[3]+x_x);
        F[1]=(1.07-x[0]-0.5*sin(x[1])-x[2]-x[3]+y_y);
        F[2]=(3.74-0.5*x[0]-x[1]-cos(x[2])-x[3]+x_x);
        F[3]=(0.79-x[0]-0.5*x[1]-x[2]-sin(x[3])+y_y);
F_1[0][0]=-0.5*sin(x[0]);
F_1[0][2]= 1;
F_1[0][3]= 1;
F_1[1][0]=1; 
F_1[1][1]=0.5*cos(x[1]);
 F_1[1][2]= 1; 
F_1[1][3]= 1;
F_1[2][0]=0.5;
F_1[2][1]=1;
F_1[2][2]=-sin(x[2]); F_1[2][3]= 1;
F_1[3][0]=1;
F_1[3][1]=0.5;
F_1[3][2]= 1;
F_1[3][3]=cos(x[3]);
DooLittle(F_1,F);
        if(MAX(x_)/MAX(x)<=1.0e-12)
            return;     
        for(j=0;j<n;j++)
        {
            x[j]+=x_[j];
        }       
    }
    cout<<"迭代2000次没有带到精度要求"<<endl;
    return;
}


void Solve_tu()//求解对应的t、u子程序
{
    int i,j;
    double x_x[11],y_y[21];   
    for(i=0;i<11;i++)
        x_x[i]=0.08*i;  
    for(j=0;j<21;j++)
        y_y[j]=0.5+0.05*j;
    for(i=0;i<11;i++)
        for(j=0;j<21;j++)
        {      
            NEWTON(x_x[i],y_y[j]);
            t[i][j]=x[0];
            u[i][j]=x[1];
        }   
    for(i=0;i<8;i++)
        x_x[i]=0.1*(i+1);   
    for(j=0;j<5;j++)
        y_y[j]=0.5+0.2*(j+1);
    for(i=0;i<8;i++)
        for(j=0;j<5;j++)
        {    
            NEWTON(x_x[i],y_y[j]);
            tt[i][j]=x[0];
            uu[i][j]=x[1];
        }
}

void Lagrange(double *t,double *u,double *U,int la,int lb,int flag)//分片二次代数插值
{
int i,j,k,m,p,c,d,q;
double a[11][21],b[11][21],Z[6][6],temp1,temp2,L1[3],L2[3];

Z[0][0]=-0.5 ; Z[0][1]=-0.34; Z[0][2]= 0.14; Z[0][3]= 0.94; Z[0][4]= 2.06;  Z[0][5]=  3.5;
Z[1][0]=-0.42; Z[1][1]=-0.5 ; Z[1][2]=-0.26; Z[1][3]= 0.3 ; Z[1][4]= 1.18;  Z[1][5]= 2.38;
Z[2][0]=-0.18; Z[2][1]=-0.5 ; Z[2][2]=-0.5 ; Z[2][3]=-0.18; Z[2][4]= 0.46;  Z[2][5]= 1.42;
Z[3][0]= 0.22; Z[3][1]=-0.34; Z[3][2]=-0.58; Z[3][3]=-0.5 ; Z[3][4]=-0.1 ;  Z[3][5]= 0.62;
Z[4][0]= 0.78; Z[4][1]=-0.02; Z[4][2]=-0.5 ; Z[4][3]=-0.66; Z[4][4]=-0.5 ;  Z[4][5]=-0.02;
Z[5][0]= 1.5 ; Z[5][1]= 0.46; Z[5][2]=-0.26; Z[5][3]=-0.66; Z[5][4]=-0.74;  Z[5][5]=-0.5; 
for(i=0;i<la;i++)
  {
      for(j=0;j<lb;j++)
      {       
        c=int(*(u+i*lb+j)/0.4);
        d=int((*(u+i*lb+j)-c*0.4)/(0.5*0.4));
        a[i][j]=(c+d)*0.4;       

        c=int(*(t+i*lb+j)/0.2);
        d=int((*(t+i*lb+j)-c*0.2)/(0.5*0.2));
        b[i][j]=(c+d)*0.2;

        if(a[i][j]<0.4)         a[i][j]=0.4;
        if(a[i][j]>1.6)         a[i][j]=1.6;
        if(b[i][j]<0.2)         b[i][j]=0.2;
        if(b[i][j]>0.8)         b[i][j]=0.8;

        for(k=0;k<3;k++)
        {   
            temp1=1;temp2=1;
            for(m=0;m<3;m++)
                if(m!=k)
                {
                    temp1*=(*(t+i*lb+j)-(b[i][j]+(m-1)*0.2))/(b[i][j]+(k-1)*0.2-(b[i][j]+(m-1)*0.2));
                    temp2*=(*(u+i*lb+j)-(a[i][j]+(m-1)*0.4))/(a[i][j]+(k-1)*0.4-(a[i][j]+(m-1)*0.4));
                }
            L1[k]=temp1;
            L2[k]=temp2;
        }

        temp1=0;
        for(k=0;k<3;k++)
            for(m=0;m<3;m++)
            {
                p=int(b[i][j]/0.2)-1+k;
                q=int(a[i][j]/0.4)-1+m;         
                temp1+=L1[k]*L2[m]*Z[p][q];
            }
        *(U+i*lb+j)=temp1;

        if(flag)
        {
         printf("%s%2d%s%.2f%s%2d%s%.2f%s%19.11e%s","x(",i,")=",0.08*i,", y(",j,")=",0.5+0.05*j,", ",*(U+i*lb+j),"; ");
         if(!fmod(j+1,2))
            cout<<endl;
         if(j==lb-1)
             cout<<endl;
        }
      }
  }
}


void Multi(double *a, double *b, double *c, int la, int lb, int lc, int r, int s, int t)

{
   int i, j, k;
   for (i=0; i<r; i++)
       for (j=0; j<t; j++)
       {
           *(c+i*lc+j)=0;
           for (k=0; k<s; k++)*(c+i*lc+j)+=*(a+i*la+k)*(*(b+k*lb+j));
       }
}

double Inverse(double *a, double *b, int la, int lb, int N)
          {
   int i, j, k;
   double temp;
   for(i=0; i<N; i++)
      for(j=0; j<N; j++)
        if (i==j)
            *(b+i*lb+j)=1;
        else
            *(b+i*lb+j)=0;
   for (k=0; k<N; k++){
      j=k;
      for (i=k+1; i<N; i++)
          if (fabs(*(a+i*la+k))>fabs(*(a+j*la+k)))    j=i;
      if (j!=k)
          for (i=0; i<N; i++)
          {
                  temp=*(a+j*la+i);
                  *(a+j*la+i)=*(a+k*la+i);
                  *(a+k*la+i)=temp;
                  temp=*(b+j*lb+i);
                  *(b+j*lb+i)=*(b+k*lb+i);
                  *(b+k*lb+i)=temp;
          }
      if (*(a+k*la+k)==0)
          return 0;
      if ((temp=*(a+k*la+k))!=1)
      for (i=0; i<N; i++){
          *(a+k*la+i)/=temp;
          *(b+k*lb+i)/=temp;
      }
      for (i=0; i<N; i++)
         if ((*(a+i*la+k)!=0) && (i!=k)){
             temp=*(a+i*la+k);
             for (j=0; j<N; j++){
                 *(a+i*la+j)-=temp*(*(a+k*la+j));
                 *(b+i*lb+j)-=temp*(*(b+k*lb+j));
             }
         }
   }
   return 0;
}

void solve_C()
{
     int    i,j,r,s;
     double t1[21][21], t2[21][21], t3[21][21],d[9][9],e[9][9];
     double B[11][9], B_T[9][11], G[21][9], G_T[9][21],P[11][21];
     double temp, FangCha;

     for(i=0;i<9;i++)
     {
           for(j=0;j<11;j++)
           {
               B[j][i]=pow(0.08*j,i);
               B_T[i][j]=pow(0.08*j,i);
           }
           for(j=0;j<21;j++)
           {
               G[j][i]=pow(0.5+0.05*j,i);
               G_T[i][j]=pow(0.5+0.05*j,i);
           }
      }   

     for (k=0; k<9; k++) 
      {
          FangCha=0;
          Multi(B_T[0], B[0], t1[0], 11, 9, 21, k+1, 11, k+1);   
          Inverse(t1[0], c[0], 21, 9, k+1);
          Multi(e[0], c[0], d[0], 9, 9, 9, k+1, k+1, k+1);
          Multi(c[0], B_T[0], t1[0], 9, 11, 21, k+1, k+1, 11);
          Multi(t1[0], U[0], t2[0], 21, 21, 21, k+1, 11, 21);
          Multi(G_T[0], G[0], t1[0], 21, 9, 21, k+1, 21, k+1);
          Inverse(t1[0], c[0], 21, 9, k+1);
          Multi(G[0], c[0], t3[0], 9, 9, 21, 21, k+1, k+1);
          Multi(t2[0], t3[0], c[0], 21, 21, 9, k+1, 21, k+1);
          for(i=0;i<11;i++)
             for(j=0;j<21;j++)
             {
                 temp=0;
                 for(r=0;r<k+1;r++)
                   for(s=0;s<k+1;s++)
                      temp+=c[r][s]*B[i][r]*G[j][s];
                  P[i][j]=temp;      
                  FangCha+=(U[i][j]-temp)*(U[i][j]-temp);
              }
           printf("%s%d%s%.11e%s","k=",k,"; Sigma=",FangCa,";\n");
          if(FangCha<=1.0e-7)
          {
             printf("%s%d%s%.11e%s"," k=",k,"; sigma=",FangCha,";\n系数c(r,s)如下：\n");
             for(i=0;i<k+1;i++)
             {
                 for(j=0;j<k+1;j++)
                 {
                     printf("%s%d%s%d%s%19.11e%s","c(",i,",",j,")=",c[i][j],"; ");
                     if(j==(k-1)/2)
                     cout<<endl;
                 }

                 cout<<endl;

             }
             cout<<endl;
          return;
        }
    }
    cout<<"经过拟合没有达到所需精度；"<<endl;
    return;
}

#define n 4
static double x[n],x_[n],U[11][21],t[11][21],u[`
1][21],tt[8][5],uu[8][5],UU[8][5],c[9][9];
int k;

void main()
{ 
        NewTown(double,double);
     DoLittle(double **,double *);
    getchar();
}