北航数值分析报告大作业第八题.doc
《北航数值分析报告大作业第八题.doc》由会员分享，可在线阅读，更多相关《北航数值分析报告大作业第八题.doc（18页珍藏版）》请在冰点文库上搜索。
实用文档
北京航空航天大学
数值分析大作业八
学院名称自动化
专业方向控制工程
学号
学生姓名许阳
教 师孙玉泉
日期2014年11月26日
一．题目
关于x,y,t,u,v,w的方程组（A.3）
（A.3）
以及关于z,t,u的二维数表（见表A-1）确定了一个二元函数z=f（x,y）。
表A-1二维数表
t
z
u
0
0.4
0.8
1.2
1.6
2
0
-0.5
-0.34
0.14
0.94
2.06
3.5
0.2
-0.42
-0.5
-0.26
0.3
1.18
2.38
0.4
-0.18
-0.5
-0.5
-0.18
0.46
1.42
0.6
0.22
-0.34
-0.58
-0.5
-0.1
0.62
0.8
0.78
-0.02
-0.5
-0.66
-0.5
-0.02
1.0
1.5
0.46
-0.26
-0.66
-0.74
-0.5
1.试用数值方法求出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.与的数表的获得
对区域上的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）上，于是选取为迭代初值。
通过以上方法求得与对应的。
（2）分片二元双二次代数插值
为保证代数插值的收敛性，应采用分片低次插值。
故此使用分片双二次代数插值。
给定如满足如下关系式：
，
，
则选择为插值节点，相应插值多项式为
其中
如果，则上式取m=1或m=4；如果或，则取n=1或n=4。
得到表达式后，直接带入，得到的值即为与对应的。
2.乘积型最小二乘曲面拟合
使用乘积型最小二乘拟合，根据k值不用，有基函数矩阵如下：
，
数表矩阵如下：
记C=[]，则系数的表达式矩阵为：
通过求解如下线性方程，即可得到系数矩阵C。
3.计算（i=1,2,…,8;j=1,2,…,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（kF（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）111;10.5*cos（u）11;0.51-sin（v）1;10.51cos（w）];
deltaX=ones（4,1）;
deltaX=dF^-1*（-1）*F;
ifmax（abs（deltaX））/abs（x）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
分片双二次代数插值子程序：
functionf=p22（t,u）
z=[-0.5-0.340.140.942.063.5;
-0.42-0.5-0.260.31.182.38;
-0.18-0.5-0.5-0.180.461.42;
0.22-0.34-0.58-0.5-0.10.62;
0.78-0.02-0.5-0.66-0.5-0.02;
1.50.46-0.26-0.66-0.74-0.5];
if（t<=0.3）
i=1;
end
if（t>0.3&&t<=0.5）
i=2;
end
if（t>0.5&&t<=0.7）
i=3;
end
if（t>0.7）
i=4;
end
ifu<=0.6
j=1;
end
ifu>0.6&&u<=1
j=2;
end
ifu>1&&u<=1.4
j=3;
end
ifu>1.4
j=4;
end
fork=i:
i+2
num=1;
fora=i:
i+2
ifa~=k
num=num*（t-0.2*（a-1））/（0.2*（k-1）-0.2*（a-1））;
end
end
l（k）=num;
end
forr=j:
j+2
num=1;
fora=j:
j+2
ifa~=r
num=num*（u-0.4*（a-1））/（0.4*（r-1）-0.4*（a-1））;
end
end
m（r）=num;
end
sum=0;
fork=i:
i+2
forr=j:
j+2
sum=sum+l（k）*m（r）*z（k,r）;
end
end
f=sum;
最小二乘曲面拟合子程序：
function[C,k,sigma]=fitxy（A）
k=1;N=10;
whilekB=ones（11,k+1）;
G=ones（21,k+1）;
fori=1:
k
forl=1:
11
B（l,i+1）=（0.08*（l-1））^i;
end
form=1:
21
G（m,i+1）=（0.5+0.05*（m-1））^i;
end
end
U=zeros（11,21）;
fori=1:
11
forj=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;
fori=1:
11
forj=1:
21
p=0;
forr=0:
k
fors=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）-p）^2;
end
End
k
sigma
ifsigma<=1e-7
break;
end
k=k+1;
end
主程序：
clc;clear;
A1=zeros（11*21,3）;
fori=1:
11
x（i）=0.08*（i-1）;
forj=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））;
end
end
[C,k,sigma]=fitxy（A1）;
A2=zeros（40,4）;
fori=1:
8
x1（i）=0.1*i;
forj=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））;
end
end
fori=1:
8
forj=1:
5
p=0;
forr=0:
k
fors=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;
end
end
A1
A2
程序结果：
数表：
k和:
数表：
标准
升级会员 