北航数值分析报告大作业第八题.doc

1、实用文档北 京 航 空 航 天 大 学数值分析大作业八学院名称 自动化 专业方向 控制工程 学 号 学生姓名 许阳 教师 孙玉泉 日 期 2014 年 11月26 日 一题目关于x , y , t , u , v , w 的方程组(A.3) (A.3)以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z=f(x , y)。表A-1 二维数表tzu00.40.81.21.620-0.5-0.340.140.942.063.50.2-0.42-0.5-0.260.31.182.380.4-0.18-0.5-0.5-0.180.461.420.60.22-0.34-0.58-0.

2、5-0.10.620.80.78-0.02-0.5-0.66-0.5-0.021.01.50.46-0.26-0.66-0.74-0.51. 试用数值方法求出f (x , y) 在区域上的近似表达式要求p(x , y)以最小的k值达到以下的精度其中。2. 计算 (i=1,2,8 ; j=1,2,5) 的值，以观察p(x , y) 逼近f (x , y)的效果，其中。二 算法设计（一）总体思路1.题目要求对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。与的数表由方程组与表A-1得到。2.与1使用相同方法求得，由计算得出的p(x,y)直接带入求得。（二）算法实现1. 与的数表的获得对区域上的

3、f (x , y)值可由方程组及二维数表得到。将区域D上的分别回代于方程组(A.3)，成为关于t,u,v,w的4元非线性方程组，解出每个对应的t,u。再通过表A-1进行插值近似，得到相应的z值。对应的z即为D区域上对应的。从而得到与的数表。(1) 4元非线性方程组求解代入(A.3)后，原方程组变为关于t,u,v,w的4元非线性方程组。观察到方程组中方程形式较为简单，易于对变量t,u,v,w求偏导数，故而选用Newton法对方程组求解。计算方程组矩阵为：计算方程组偏导数矩阵为：迭代公式为：，k=0,1,2,n其中为线性方程组的解。取为迭代终止条件。由表A-1观察到t,u基本在（0,2）上，于是选

4、取为迭代初值。通过以上方法求得与对应的。(2) 分片二元双二次代数插值为保证代数插值的收敛性，应采用分片低次插值。故此使用分片双二次代数插值。给定如满足如下关系式：，则选择为插值节点，相应插值多项式为其中 如果，则上式取m=1或m=4；如果或，则取n=1或n=4。得到表达式后，直接带入，得到的值即为与对应的。2. 乘积型最小二乘曲面拟合使用乘积型最小二乘拟合，根据k值不用，有基函数矩阵如下： ， 数表矩阵如下：记C=，则系数的表达式矩阵为：通过求解如下线性方程，即可得到系数矩阵C。3. 计算 (i=1,2,8 ; j=1,2,5) 的值的计算与相同。将代入原方程组，求解响应进行分片双二次插值求

5、得。的计算则可以直接将代入所求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(kN) 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)

6、1;1 0.5 1 cos(w); deltaX=ones(4,1); deltaX=dF-1*(-1)*F; if max(abs(deltaX)/abs(x)ep t=t+deltaX(1,1); u=u+deltaX(2,1); v=v+deltaX(3,1); w=w+deltaX(4,1); break; end t=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

7、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 (t0.3&t0.5&t0.7) i=4;end if u0.6&u1&u1.4 j=4;end for k=i:i+2 num=1; for a=i:i+2 if a=k num=num*(t-0.2*(a-1)/(0.2*(k-1)-0.2*(a-1)

8、; end end l(k)=num; endfor r=j:j+2 num=1; for a=j:j+2 if a=r num=num*(u-0.4*(a-1)/(0.4*(r-1)-0.4*(a-1); end end m(r)=num; endsum=0;for k=i:i+2 for r=j:j+2 sum=sum+l(k)*m(r)*z(k,r); endendf=sum;最小二乘曲面拟合子程序：function C,k,sigma=fitxy(A)k=1;N=10;while kN B=ones(11,k+1); G=ones(21,k+1); for i=1:k for l=1:

9、11 B(l,i+1)=(0.08*(l-1)i; end for m=1:21 G(m,i+1)=(0.5+0.05*(m-1)i; end end U=zeros(11,21); for i=1:11 for j=1:21 U(i,j)=A(i-1)*21+j,3); end end C=(B*B)-1*B*U*G*(G*G)-1; sigma=0; for i=1:11 for j=1:21 p=0; for r=0:k for s=0:k p=p+C(r+1,s+1)*(0.08*(i-1)r*(0.5+(0.05*(j-1)s; end end sigma=sigma+(U(i,j)

10、-p)2; endEndksigma if sigma=1e-7 break; end k=k+1;end主程序：clc;clear;A1=zeros(11*21,3);for i=1:11 x(i)=0.08*(i-1); for j=1:21 y(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

11、 i=1:8 x1(i)=0.1*i; for j=1:5 y1(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); endend for i=1:8 for j=1:5 p=0; for r=0:k for s=0:k p=p+C(r+1,s+1)*(0.1*i)r*(0.5+(0.2*j)s; end end A2(i-1)*5+j,4)=p; endendA1A2程序结果：数表：k和:数表：标准