北京航空航天大学数值分析大作业八学院名称自动化专业方向控制工程学号学生姓名许阳教师孙玉泉日期2014 年11月26 日一．题目关于x , y , t , u , v , w 的方程组(A.3)⎪⎪⎩⎪⎪⎨⎧=-+++=-+++=-+++=-+++79.0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。

表A-1 二维数表1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式∑∑===k i kj s r rs y x c y x p 00),(要求p (x , y )以最小的k 值达到以下的精度∑∑==-≤-=100207210)],(),([i j i i i i y x p y x f σ其中j y i x i i 05.05.0,08.0+==。

2. 计算),(),,(****j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。

二．算法设计（一）总体思路1.题目要求∑∑===ki kj s r rs y x c y x p 00),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。

),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。

2.),(**j i y x f 与1使用相同方法求得，),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。

（二）算法实现1. ),(i i y x 与),(i i y x f 的数表的获得对区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的f (x , y )值可由方程组及二维数表得到。

将区域D 上的),(i i y x 分别回代于方程组(A.3)，成为关于t,u,v,w 的4元非线性方程组，解出每个),(j i y x 对应的t,u 。

再通过表A-1进行插值近似，得到相应的z 值。

对应的z 即为D 区域上),(j i y x 对应的),(j i y x f 。

从而得到),(j i y x 与),(i i y x f 的数表。

(1) 4元非线性方程组求解),(i i y x 代入(A.3)后，原方程组变为关于t,u,v,w 的4元非线性方程组。

观察到方程组中方程形式较为简单，易于对变量t,u,v,w 求偏导数，故而选用Newton 法对方程组求解。

计算方程组矩阵为：⎪⎪⎪⎪⎪⎭⎫⎝⎛--+++--+++--+++--+++=79.0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0),,,(i i i i y w v u t x w v u t y w v u t x w v u t w v u t F计算方程组偏导数矩阵为：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=w v u tw v u t F cos 15.011sin 15.011cos 5.01111sin 5.0),,,('迭代公式为：⎪⎪⎪⎪⎪⎭⎫⎝⎛∆∆∆∆+⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛++++)()()()()()()()()1()1()1()1(k k k k k k k k k k k k w v u t w v u t w v u t ，k=0,1,2,…,n 其中⎪⎪⎪⎪⎪⎭⎫ ⎝⎛∆∆∆∆)()()()(k k k k w v u t 为线性方程组),,,(),,,(')()()()()()()()()()()()(k k k k k k k k k k k k w v u t F w v u t w v u t F -=⎪⎪⎪⎪⎪⎭⎫⎝⎛∆∆∆∆的解。

取ε≤∆∆∆∆|)||,||,||,m ax (|/|)||,||,||,m ax (|)()()()()()()()(k k k k k k k k w v u t w v u t 为迭代终止条件。

由表A-1观察到t,u 基本在（0,2）上，于是选取)1,1,1,1(),,,)0()0()0()0(=w v u t （为迭代初值。

通过以上方法求得与),(j i y x 对应的),t (ij ij u 。

(2) 分片二元双二次代数插值为保证代数插值的收敛性，应采用分片低次插值。

故此使用分片双二次代数插值。

)5,...,1,0;5,...,1,0(4.0,2.0====n m n u m t n m给定),t (ij ij u 如满足如下关系式：1.01.0+≤<-m ij m t t t ，41<<m2.02.0+≤<-n ij n u u u ，41<<n则选择)2,1,;2,1,)(,(++=++=n n n r m m m k u t r k 为插值节点，相应插值多项式为∑∑+=+==2222),()()(),(m m k n n r r k ij r ij k ij ij u t z u l t l u t p其中∏+≠=--=2)(m k q m q qkq ij ij k t t t t t l )21,++=m m m k ，（∏+≠=--=2)(n kq n q qk q ij ij r u u u u u l )2,1,++=n n n r （如果7.03.0>≤ij ij t t 或，则上式取m=1或m=4；如果6.0≤ij u 或4.1>ij u ，则取n=1或n=4。

得到),(22ij ij u t p 表达式后，直接带入),t (ij ij u ，得到的值即为与),(j i y x 对应的),(j i y x f 。

2. 乘积型最小二乘曲面拟合使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下：⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ]，则系数rs c 的表达式矩阵为：11-)(-=G G UG B B B C T TT ）（通过求解如下线性方程，即可得到系数矩阵C 。

UG B G G C B B T T T =)()(3.计算),(),,(****j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值),(**j i y x f 的计算与),(j i y x f 相同。

将),(**j i y x 代入原方程组，求解响应),(**ij ij u t 进行分片双二次插值求得),(**j i y x f 。

),(**j i y x p 的计算则可以直接将),(**j i y x 代入所求p(x,y)。

三．Matlab源程序及结果牛顿法解非线性方程组子程序：function [t,u]=newt(x,y)t=1; u=1; v=1; w=1;ep=1e-12;k=1;N=100;while(k<N)F(1,1)=0.5*cos(t)+u+v+w-x-2.67;F(2,1)=t+0.5*sin(u)+v+w-y-1.07;F(3,1)=0.5*t+u+cos(v)+w-x-3.74;F(4,1)=t+0.5*u+v+sin(w)-y-0.79;dF=[-0.5*sin(t) 1 1 1;1 0.5*cos(u) 1 1;0.5 1 -sin(v) 1;1 0.5 1 cos(w)]; deltaX=ones(4,1);deltaX=dF^-1*(-1)*F;if max(abs(deltaX))/abs(x)<ept=t+deltaX(1,1);u=u+deltaX(2,1);v=v+deltaX(3,1);w=w+deltaX(4,1);break;endt=t+deltaX(1,1);u=u+deltaX(2,1);v=v+deltaX(3,1);w=w+deltaX(4,1);k=k+1;end分片双二次代数插值子程序：function f=p22(t,u)z=[-0.5 -0.34 0.14 0.94 2.06 3.5;-0.42 -0.5 -0.26 0.3 1.18 2.38;-0.18 -0.5 -0.5 -0.18 0.46 1.42;0.22 -0.34 -0.58 -0.5 -0.1 0.62;0.78 -0.02 -0.5 -0.66 -0.5 -0.02;1.5 0.46 -0.26 -0.66 -0.74 -0.5];if (t<=0.3)i=1;endif (t>0.3&&t<=0.5)i=2;endif (t>0.5&&t<=0.7)i=3;endif (t>0.7)i=4;endif u<=0.6j=1;endif u>0.6&&u<=1j=2;endif u>1&&u<=1.4j=3;endif u>1.4j=4;endfor k=i:i+2num=1;for a=i:i+2if a~=knum=num*(t-0.2*(a-1))/(0.2*(k-1)-0.2*(a-1));endendl(k)=num;endfor r=j:j+2num=1;for a=j:j+2if a~=rnum=num*(u-0.4*(a-1))/(0.4*(r-1)-0.4*(a-1));endendm(r)=num;endsum=0;for k=i:i+2for r=j:j+2sum=sum+l(k)*m(r)*z(k,r);endendf=sum;最小二乘曲面拟合子程序：function [C,k,sigma]=fitxy(A)k=1;N=10;while k<NB=ones(11,k+1);G=ones(21,k+1);for i=1:kfor l=1:11B(l,i+1)=(0.08*(l-1))^i;endfor m=1:21G(m,i+1)=(0.5+0.05*(m-1))^i;endendU=zeros(11,21);for i=1:11for j=1:21U(i,j)=A((i-1)*21+j,3);endendC=(B'*B)^-1*B'*U*G*(G'*G)^-1;sigma=0;for i=1:11for j=1:21p=0;for r=0:kfor s=0:kp=p+C(r+1,s+1)*(0.08*(i-1))^r*(0.5+(0.05*(j-1)))^s;endendsigma=sigma+(U(i,j)-p)^2;endEndksigmaif sigma<=1e-7break;endk=k+1;end主程序：clc;clear;A1=zeros(11*21,3);for i=1:11x(i)=0.08*(i-1);for j=1:21y(j)=0.5+0.05*(j-1);[t(i,j),u(i,j)]=newt(x(i),y(j));A1((i-1)*21+j,1)=x(i);A1((i-1)*21+j,2)=y(j);A1((i-1)*21+j,3)=p22(t(i,j),u(i,j));endend[C,k,sigma]=fitxy(A1);A2=zeros(40,4);for i=1:8x1(i)=0.1*i;for j=1:5y1(j)=0.5+0.2*j;[t1(i,j),u1(i,j)]=newt(x1(i),y1(j));A2((i-1)*5+j,1)=x1(i);A2((i-1)*5+j,2)=y1(j);A2((i-1)*5+j,3)=p22(t1(i,j),u1(i,j));endendfor i=1:8for j=1:5p=0;for r=0:kfor s=0:kp=p+C(r+1,s+1)*(0.1*i)^r*(0.5+(0.2*j))^s;endendA2((i-1)*5+j,4)=p;endendA1A2程序结果：数表：)),(,,i i i i y x f y x （k和 :数表)),(),,(,,(******j i j i j i y x p y x f y x ：。