一、题目
关于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)=∑s=0kxrys
要求p(x,y)一最小的k值达到以下的精度
σ=∑i=010∑j=020[f(xi−yj)−p(xi,yj)]2≤10−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)
代入非线性方程组中,用牛顿法解出
ti
和
uj
;
2.以采取分片二次插值,选择(m,n)满足
tm−h/2<ti≤tm+h/2,2≤m≤3
un−τ/2<uj≤un+τ/2,2≤n≤3
如果
ti≤t1+h/2
或
ti>t4−h/2
,则m=1或4;如果
uj≤u1+τ/2
或
uj>u4−τ/2
,则n=1或n=4。选择
(tk,ur)(k=m−1,m,m+1;r=n−1,n,n+1)
为插值节点,相应的Lagrange形式的插值多项式为
p22(t,u)=∑k=m−1m+1∑r=n−1n+1lk(t)l̂ rf(tk,ur)
其中
lk(t)=∏w=m−1;w≠rm+1t−twtk−tw(k=m−1,m,m+1)
lr(u)=∏w=n−1;w≠rn+1y−ywyr−yw(r=n−1,n,n+1)
并将 tI 和 uj 代入 p22(t,u) ,便得到了数表 xi,yj,f(xi,yj) 。
3.进行曲面拟和系数矩阵 C=[crs∗] , C=(BTB)−1BTUG(GTG)−1
其中
B=[ϕc(xi)]=⎡⎣⎢⎢⎢⎢11M1x0x1Mx10LLMLxkkxk1Mxk10⎤⎦⎥⎥⎥⎥
,
G=[ψs(yj)]=⎡⎣⎢⎢⎢⎢11M1y0y1My10LLMLykkyk1Myk10⎤⎦⎥⎥⎥⎥
U=[f(xi,yj)]
k从0逐渐增大,直到 σ≤10−7 ,便得到了要求精度的系数 crs 。
三、运算结果
得到插值节点:
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();
}