北 京 航 空 航 天 大 学
数值分析大作业八
学院名称 自动化 专业方向 控制工程 学 号 学生姓名 许阳 教
师 孙玉泉
日 期 2014 年 11月26 日
一．题目
关于x , y , t , u , v , w 的方程组(A.3)
?0.5cost?u?v?w?x?2.67?t?0.5sinu?cosv?w?y?1.07? (A.3) ?0.5t?u?cosv?w?x?3.74???t?0.5u?v?sinw?y?0.79以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z=f(x , y)。
表A-1 二维数表
t z u 0 0.2 0.4 0.6 0.8 1.0 0 0.4 0.8 1.2 1.6 2 -0.5 -0.42 -0.18 0.22 0.78 1.5 -0.34 -0.5 -0.5 -0.34 -0.02 0.46 0.14 -0.26 -0.5 -0.58 -0.5 -0.26
0.94 0.3 -0.18 -0.5 -0.66 -0.66 2.06 1.18 0.46 -0.1 -0.5 -0.74 3.5 2.38 1.42 0.62 -0.02 -0.5 （ x,y)|0?x?0.8,0.5?y?1.5}上的1. 试用数值方法求出f (x , y) 在区域D?{近似表达式
p(x,y)???crsxrys
i?0j?0kk要求p(x , y)以最小的k值达到以下的精度
????[f(xi,yi)?p(xi,yi)]2?10?7
i?0j?01020其中xi?0.08i,yi?0.5?0.05j。
**2. 计算f(xi*,y*j),p(xi,yj) (i=1,2,…,8 ; j=1,2,…,5) 的值，以观察p(x , y) 逼
近f (x , y)的效果，其中xi*?0.1i,y*j?0.5?0.2j。
二．算法设计
（一）总体思路
1.题目要求p(x,y)???crsxrys对f(x, y) 进行拟合，可选用乘积型最小二乘
i?0j?0kk拟合。(xi,yi)与f(xi,yi)的数表由方程组与表A-1得到。
**2.f(xi*,y*j)与1使用相同方法求得，p(xi,yj)由计算得出的p(x,y)直接带入
(xi*,y*j)求得。
（二）算法实现
1. (xi,yi)与f(xi,yi)的数表的获得
（ x,y)|0?x?0.8,0.5?y?1.5}上的f (x , y)值可由方程组及二维对区域D?{数表得到。将区域D上的(xi,yi)分别回代于方程组(A.3)，成为关于t,u,v,w的4元非线性方程组，解出每个(xi,yj)对应的t,u。再通过表A-1进行插值近似，得到相应的z值。对应的z即为D区域上(xi,yj)对应的f(xi,yj)。从而得到(xi,yj)与f(xi,yi)的数表。
(1) 4元非线性方程组求解
(xi,yi)代入(A.3)后，原方程组变为关于t,u,v,w的4元非线性方程组。观察到方程组中方程形式较为简单，易于对变量t,u,v,w求偏导数，故而选用Newton法对方程组求解。
计算方程组矩阵为：
?0.5cost?u?v?w?xi?2.67????t?0.5sinu?cosv?w?yi?1.07?F(t,u,v,w)??
0.5t?u?cosv?w?xi?3.74????t?0.5u?v?sinw?y?0.79?i??
计算方程组偏导数矩阵为：
111???0.5sint??10.5cosu11?? F'(t,u,v,w)??0.51?sinv1????10.51cosw???
迭代公式为：
?t(k?1)??t(k)???t(k)??(k?1)??(k)??(k)??u??u???u??v(k?1)???v(k)????v(k)?，k=0,1,2,…,n ???????w(k?1)??w(k)???w(k)???????其中
??t(k)???t(k)??(k)??(k)????u?(k)(k)(k)(k)??u(k)(k)(k)(k)F'(t,u,v,w)??F(t,u,v,w)的解。 为线性方程组??v(k)??(k)???v?????w(k)???w(k)?????取max(|?t(k)|,|?u(k)|,|?v(k)|,|?w(k)|)/max(|t(k)|,|u(k)|,|v(k)|,|w(k)|)??为迭代终止条件。
由表A-1观察到t,u基本在（0,2）上，于是选取（t(0),u(0),v(0),w(0))?(1,1,1,1)为迭代初值。
通过以上方法求得与(xi,yj)对应的(tij,uij)。 (2) 分片二元双二次代数插值
为保证代数插值的收敛性，应采用分片低次插值。故此使用分片双二次代数插值。
tm?0.2m,un?0.4n(m?0,1,...,5;n?0,1,...,5)
给定(tij,uij)如满足如下关系式：
tm?0.1?tij?tm?0.1，1?m?4
un?0.2?uij?un?0.2，1?n?4
则选择(tk,ur)(k?m,m?1,m?2;r?n,n?1,n?2)为插值节点，相应插值多项式为
p22(tij,uij)???lk(tij)lr(uij)z(tk,ur)
k?mr?nm?2n?2其中
lk(tij)??m?2tij?tqt?tq （k?m,m?1，m?2)
q?mkq?kn?2lr(uij)??q?nq?kuij?uquk?uq （r?n,n?1,n?2)
如果tij?0.3或tij?0.7，则上式取m=1或m=4；如果uij?0.6或uij?1.4，则取n=1或n=4。
得到p22(tij,uij)表达式后，直接带入(tij,uij)，得到的值即为与(xi,yj)对应的
f(xi,yj)。
2. 乘积型最小二乘曲面拟合
使用乘积型最小二乘拟合，根据k值不用，有基函数矩阵如下：
0k0k?x0??y0??x0?y0????B?????? ， G??????
?x0?xk??y0?yk?iij????j数表矩阵如下：
?f(x0,y0)??U?????f(x,y)?i0?f(x0,yj)???? f(xi,yj)??记C=[crs]，则系数crs的表达式矩阵为：
-1TC?（BTB）BUG(GTG)?1
通过求解如下线性方程，即可得到系数矩阵C。
(BTB)C(GTG)?BTUG