北航数值分析报告大作业第八题.pdf
实用文档北北 京京 航航 空空 航航 天天 大大 学学数值分析大作业八数值分析大作业八学院名称学院名称自动化自动化专业方向专业方向控制工程控制工程学学号号学生姓名学生姓名许阳许阳教教师师孙玉泉孙玉泉日日期期 20142014 年年 1111 月月 2626 日日标准实用文档一题目一题目关于 x,y,t,u,v,w 的方程组(A.3)0.5cost u v w x 2.67t 0.5sinu cosv w y 1.07(A.3)0.5t u cosv w x 3.74t 0.5u vsinw y 0.79以及关于 z,t,u 的二维数表（见表 A-1）确定了一个二元函数 z=f(x,y)。表表 A-1A-1 二维数表二维数表tzu00.20.40.60.81.000.40.81.21.62-0.5-0.42-0.180.220.781.5-0.34-0.5-0.5-0.34-0.020.460.14-0.26-0.5-0.58-0.5-0.260.940.3-0.18-0.5-0.66-0.662.061.180.46-0.1-0.5-0.743.52.381.420.62-0.02-0.51.试用数值方法求出f(x,y)在区域D （x,y)|0 x 0.8,0.5 y 1.5上的近似表达式p(x,y)crsxrysi0 j0kk要求 p(x,y)以最小的 k 值达到以下的精度 f(xi,yi)p(xi,yi)2107i0 j01020其中xi 0.08i,yi 0.50.05 j。*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*j0.50.2j。标准实用文档二算法设计二算法设计（一）总体思路（一）总体思路1.题目要求p(x,y)crsxrys对 f(x,y)进行拟合，可选用乘积型最小二乘i0 j0kk拟合。(xi,yi)与f(xi,yi)的数表由方程组与表 A-1 得到。*2.f(xi*,y*j)与 1 使用相同方法求得，p(xi,yj)由计算得出的 p(x,y)直接带入(xi*,y*j)求得。（二）算法实现（二）算法实现1.1.(xi,yi)与与f(xi,yi)的数表的获得的数表的获得对区域D （x,y)|0 x 0.8,0.5 y 1.5上的 f(x,y)值可由方程组及二维数表得到。将区域 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(1)4 元非线性方程组求解元非线性方程组求解(xi,yi)代入(A.3)后，原方程组变为关于 t,u,v,w的 4 元非线性方程组。观察到方程组中方程形式较为简单，易于对变量t,u,v,w求偏导数，故而选用Newton 法对方程组求解。计算方程组矩阵为：0.5cost uvw xi2.67t 0.5sinucosvw yi1.07F(t,u,v,w)0.5t ucosvw xi3.74t 0.5uvsinw y 0.79i标准实用文档计算方程组偏导数矩阵为：1110.5sint10.5cosu11F(t,u,v,w)0.51sinv110.51cosw迭代公式为：t(k1)t(k)t(k)(k1)(k)(k)uuuv(k1)v(k)v(k)，k=0,1,2,nw(k1)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)的解。(k)(k)vvw(k)w(k)取max(|t(k)|,|u(k)|,|v(k)|,|w(k)|)/max(|t(k)|,|u(k)|,|v(k)|,|w(k)|)为迭代终止条件。（t(0),u(0),v(0),w(0)(1,1,1,1)为由表 A-1 观察到 t,u 基本在（0,2）上，于是选取迭代初值。通过以上方法求得与(xi,yj)对应的(tij,uij)。(2)(2)分片二元双二次代数插值分片二元双二次代数插值为保证代数插值的收敛性，应采用分片低次插值。故此使用分片双二次代数插值。tm 0.2m,un 0.4n(m 0,1,.,5;n 0,1,.,5)给定(tij,uij)如满足如下关系式：tm0.1 tij tm0.1，1 m 4标准实用文档un0.2 uij un0.2，1 n 4则选择(tk,ur)(k m,m1,m 2;r n,n1,n 2)为插值节点，相应插值多项式为p22(tij,uij)lk(tij)lr(uij)z(tk,ur)km rnm2 n2其中lk(tij)m2tijtqt tq（k m,m1，m2)qmkqkn2lr(uij)qnqkuijuqukuq（r n,n1,n2)如果tij 0.3或tij 0.7，则上式取 m=1 或 m=4；如果uij 0.6或uij1.4，则取 n=1 或 n=4。得到p22(tij,uij)表达式后，直接带入(tij,uij)，得到的值即为与(xi,yj)对应的f(xi,yj)。2.2.乘积型最小二乘曲面拟合乘积型最小二乘曲面拟合使用乘积型最小二乘拟合，根据 k 值不用，有基函数矩阵如下：0k0kx0 y0 x0y0B ，G x0 xky0ykiijj数表矩阵如下：f(x0,y0)U f(x,y)i0f(x0,yj)f(xi,yj)记 C=crs，则系数crs的表达式矩阵为：-1TC （BTB）B UG(GTG)1通过求解如下线性方程，即可得到系数矩阵 C。(BTB)C(GTG)BTUG标准实用文档*3.3.计算计算f(xi*,y*j),p(xi,yj)(i i=1,2,=1,2,8;,8;j j=1,2,=1,2,5),5)的值的值*f(xi*,y*j)的计算与f(xi,yj)相同。将(xi,yj)代入原方程组，求解响应(tij,uij)*进行分片双二次插值求得f(xi*,y*p(xi*,y*j)。j)的计算则可以直接将(xi,yj)代入所求 p(x,y)。标准实用文档三三MatlabMatlab 源程序及结果源程序及结果牛顿法解非线性方程组子程序：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)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 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;endif u0.6&u1&u1.4 j=4;endfor k=i:i+2 num=1;for a=i:i+2if a=k num=num*(t-0.2*(a-1)/(0.2*(k-1)-0.2*(a-1);endend l(k)=num;endfor r=j:j+2 num=1;for a=j:j+2if a=r num=num*(u-0.4*(a-1)/(0.4*(r-1)-0.4*(a-1);endend m(r)=num;endsum=0;for k=i:i+2for r=j:j+2 sum=sum+l(k)*m(r)*z(k,r);endend标准实用文档f=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:kfor l=1:11 B(l,i+1)=(0.08*(l-1)i;endfor m=1:21 G(m,i+1)=(0.5+0.05*(m-1)i;endend U=zeros(11,21);for i=1:11for j=1:21 U(i,j)=A(i-1)*21+j,3);endend C=(B*B)-1*B*U*G*(G*G)-1;sigma=0;for i=1:11for j=1:21 p=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;endend sigma=sigma+(U(i,j)-p)2;endEndksigmaif sigma=1e-7break;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);endendC,k,sigma=fitxy(A1);A2=zeros(40,4);for 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);endendfor i=1:8for j=1:5 p=0;for r=0:kfor s=0:k p=p+C(r+1,s+1)*(0.1*i)r*(0.5+(0.2*j)s;endend A2(i-1)*5+j,4)=p;endendA1A2标准实用文档程序结果：（xi,yi,f(xi,yi)数表：标准实用文档标准实用文档标准实用文档标准实用文档标准实用文档标准实用文档标准实用文档k和:标准实用文档*数表(xi*,y*j,f(xi,yj),p(xi,yj)：标准
