北航数值分析报告大作业第八题Word下载.doc
- 文档编号:13036283
- 上传时间:2022-10-03
- 格式:DOC
- 页数:18
- 大小:426.89KB
北航数值分析报告大作业第八题Word下载.doc
《北航数值分析报告大作业第八题Word下载.doc》由会员分享,可在线阅读,更多相关《北航数值分析报告大作业第八题Word下载.doc(18页珍藏版)》请在冰豆网上搜索。
0.94
2.06
3.5
0.2
-0.42
-0.26
0.3
1.18
2.38
-0.18
0.46
1.42
0.6
0.22
-0.58
-0.1
0.62
0.78
-0.02
-0.66
1.0
1.5
-0.74
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(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)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)<
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
分片双二次代数插值子程序:
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;
if(t>
0.3&
&
t<
=0.5)
i=2;
0.5&
=0.7)
i=3;
0.7)
i=4;
ifu<
=0.6
j=1;
ifu>
0.6&
u<
=1
j=2;
1&
=1.4
j=3;
1.4
j=4;
fork=i:
i+2
num=1;
fora=i:
ifa~=k
num=num*(t-0.2*(a-1))/(0.2*(k-1)-0.2*(a-1));
end
l(k)=num;
forr=j:
j+2
fora=j:
ifa~=r
num=num*(u-0.4*(a-1))/(0.4*(r-1)-0.4*(a-1));
m(r)=num;
sum=0;
forr=j:
sum=sum+l(k)*m(r)*z(k,r);
f=sum;
最小二乘曲面拟合子程序:
function[C,k,sigma]=fitxy(A)
N=10;
whilek<
N
B=ones(11,k+1);
G=ones(21,k+1);
fori=1:
k
forl=1:
11
B(l,i+1)=(0.08*(l-1))^i;
form=1:
21
G(m,i+1)=(0.5+0.05*(m-1))^i;
U=zeros(11,21);
forj=1:
U(i,j)=A((i-1)*21+j,3);
C=(B'
*B)^-1*B'
*U*G*(G'
*G)^-1;
sigma=0;
p=0;
forr=0:
fors=0:
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
sigma
ifsigma<
=1e-7
主程序:
clc;
clear;
A1=zeros(11*21,3);
fori=1:
x(i)=0.08*(i-1);
forj=1:
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));
[C,k,sigma]=fitxy(A1);
A2=zeros(40,4);
8
x1(i)=0.1*i;
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));
p=0;
forr=0:
fors=0:
p=p+C(r+1,s+1)*(0.1*i)^r*(0.5+(0.2*j))^s;
A2((i-1)*5+j,4)=p;
A1
A2
程序结果:
数表:
k和:
标准
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 报告 作业 第八
![提示](https://static.bdocx.com/images/bang_tan.gif)