北航数值分析大作业3.docx
《北航数值分析大作业3.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业3.docx（30页珍藏版）》请在冰豆网上搜索。
北航数值分析大作业3
《数值分析》计算实习题目三
一、算法设计方案
观察可知：
要求利用给定的数表进行二元函数的拟合。
由于没有直接所求拟合函数的自变量（x,y）同其函数值z的关系，而是利用数表给定另外两个变量（u,t）和z的关系，同时通过一个非线性方程组给出了（x,y）和（u,t）的对应关系。
所以要进行函数拟合之前，应先利用给定的数表和非线性方程组插值出一个关于（x,y,z）的矩形列表，然后再将该数表看作准确值进行拟合。
主要算法有：
非线性方程组的求解、二元分片二次插值和二元拟合。
其中还包括：
Gauss列主元消去法求解线性方程组，求正交多项式系等算法。
现分别列于下：
1.非线性方程组的求解
利用牛顿法求解非线性方程组。
将（x,y）作为已知量，求解关于参数（t,u,v,w）的线性方程组。
程序中将（t,u,v,w）定义为存放在X[4]数组中，再将f（t,u,v,w）=B方程转化为F（x[0],x[1],x[2],x[3]）=0,并对各自变量X[
]（
=1～4）求导，并得到关于
的线性方程组：
，通过
迭代得到非线性方程组的解。
在求解过程中设置求解精度
。
具体算法如下：
（1）给定初值
给定精度水平eps=1.0e-12
（2）利用
值得到关于
的线性方程组的系数矩阵，并求解该线性方程组，（k=0,1,2……1000）；
（3）实现迭代
，（k=0,1,2……100）；
（4）求迭代精度
，并判断EPS是否小于求解精度eps，如大于该精度返回
（2）执行下一次迭代，如小于求解精度，则程序收敛返回x1和x2的值，即所求t和u的值。
2.Gauss列主元消去法求线性方程组
（1）消元过程
对于
执
a）选行号
，使
。
b）交换
c）对于
计算
（2）回代过程
3．二元函数分片二次插值
通过求解非线性方程组后，对应每一对（x,y）可得到一个（u,t）,若要得到对应的z，需根据原有（（t,u）,z）的6×6的矩形数表和求出的（u,t）利用二元分片二次插值求出z。
对于等距节点，
，
已知。
对给定的
，若
，则选
；若
，则选
；若
，则选
。
对于等距节点，
，等上述方法一样。
其中：
4.二元函数拟合
已通过求解非线性方程组和插值得到的11×21的（（x,y）,z）二元矩形数表为基准，进行给定不同次数的二元函数最小二乘法拟合。
拟合过程选取正交函数系作为基函数，即拟合函数可表示为：
，其中
和
分别为次数为r和s的正交多项式系。
将其写为矩阵形式为：
其中：
同理：
所以：
其中
。
利用x和y的正交函数系同其幂函数系数阵以及
（2）中所求系数矩阵，最终得到所求系数矩阵，即：
。
该子程序中的矩阵乘积运算调用了通用的矩阵相乘的子程序。
二、全部源程序
#include"stdio.h"
#include"math.h"
intp1,q1;//全局变量
//二元函数插值
//x[]为插值表的x坐标点，y[]为插值表的y坐标点，z[]为插值表的函数值
//u为待插值点的x坐标，v为为待插值点的y坐标
doublecz（doublex[],doubley[],doublez[],doubleu,doublev）
{
doubleh,f,hh,w;
doublebb[10];
inti,j,ii,jj,k,k1,k2,kk;
h=0.4;f=0.2;hh=0.0;
if（u<=x[1]+h/2）
k1=1;
elseif（u>x[4]-h/2）k1=4;
else
{
for（i=2;i<=3;i++）
if（u>（x[i]-h/2）&&u<=（x[i]+h/2））
k1=i;
}
if（v<=y[1]+f/2）
k2=1;
elseif（v>y[4]-f/2）
k2=4;
else
{
for（j=2;j<=3;j++）
if（v>（y[j]-f/2）&&v<=（y[j]+f/2））
k2=j;
}
for（ii=k1-1;ii{
bb[ii]=0.0;
for（jj=k2-1;jj{
k=ii*6+jj;hh=z[k];
for（kk=k2-1;kkif（kk!
=jj）
hh*=（v-y[kk]）/（y[jj]-y[kk]）;
bb[ii]=bb[ii]+hh;
}
}
w=0.0;
for（ii=k1-1;ii<=k1+1;ii++）
{
hh=bb[ii];
for（jj=k1-1;jj<=k1+1;jj++）
if（jj!
=ii）
hh*=（u-x[jj]）/（x[ii]-x[jj]）;
w+=hh;
}
return（w）;
}
//非线性方程组牛顿迭代的Jacobi矩阵
//a表示Jacobi矩阵，X为解向量
dnewta1（X,a）
doubleX[4],a[4][4];
{
a[0][0]=-0.5*sin（X[0]）;a[0][1]=1;a[0][2]=1;a[0][3]=1;
a[1][0]=1;a[1][1]=0.5*cos（X[1]）;a[1][2]=1;a[1][3]=1;
a[2][0]=0.5;a[2][1]=1;a[2][2]=-sin（X[2]）;a[2][3]=1;
a[3][0]=1;a[3][1]=0.5;a[3][2]=1;a[3][3]=cos（X[3]）;
return;
}
//非线性方程组牛顿迭代的值向量
//X为解向量，b为值向量，xx为
dnewta2（X,b,xx,yy）
doubleX[],b[],xx[],yy[];
{
b[0]=-0.5*cos（X[0]）-X[1]-X[2]-X[3]+xx[p1]+2.67;
b[1]=-X[0]-0.5*sin（X[1]）-X[2]-X[3]+yy[q1]+1.07;
b[2]=-0.5*X[0]-X[1]-cos（X[2]）-X[3]+xx[p1]+3.74;
b[3]=-X[0]-0.5*X[1]-X[2]-sin（X[3]）+yy[q1]+0.79;
return;
}
//求向量的无穷范数
//array为向量地址，n为向量维数
doublecond（doublearray[],intn）
{
inti;
doublefan=1e-32;
for（i=0;i{
if（fabs（array[i]）>fan）
{
fan=fabs（array[i]）;
}
}
returnfan;
}//矩阵转置
//a为待转置矩阵，维数为m*n，转置后存放在矩阵b中
zhuanzhi（doublea[],doubleb[],intm,intn）
{
inti,j;
for（i=0;i{
for（j=0;j{
b[m*j+i]=a[i*n+j];
}
}
return;
}//矩阵相乘
//a矩阵维数为m*n，b矩阵维数为n*k，乘积结果存放在c矩阵中
voidbrmul（a,b,m,n,k,c）
intm,n,k;
doublea[],b[],c[];
{inti,j,l,u;
for（i=0;i<=m-1;i++）
for（j=0;j<=k-1;j++）
{u=i*k+j;c[u]=0.0;
for（l=0;l<=n-1;l++）
c[u]=c[u]+a[i*n+l]*b[l*k+j];
}
return;
}
//对矩阵部分求逆
//a为待求逆矩阵，维数为n1*n1，对其中n*n进行操作
intbrinv（doubleA[],intn,intn1）
{
double*a,d,p;
int*is,*js,i,j,k,l,u,v;
a=malloc（n*n*sizeof（double））;
is=malloc（n*sizeof（int））;
js=malloc（n*sizeof（int））;
for（i=0;i<=n-1;i++）
for（j=0;j<=n-1;j++）
*（a+i*n+j）=*（A+i*n1+j）;
for（k=0;k<=n-1;k++）
{d=0.0;
for（i=k;i<=n-1;i++）
for（j=k;j<=n-1;j++）
{l=i*n+j;p=fabs（a[l]）;
if（p>d）{d=p;is[k]=i;js[k]=j;}
}
if（d+1.0==1.0）
{free（is）;free（js）;
free（a）;
printf（"err**notinv\n"）;
return（0）;
}
if（is[k]!
=k）
for（j=0;j<=n-1;j++）
{u=k*n+j;v=is[k]*n+j;
p=a[u];a[u]=a[v];a[v]=p;
}
if（js[k]!
=k）
for（i=0;i<=n-1;i++）
{u=i*n+k;v=i*n+js[k];
p=a[u];a[u]=a[v];a[v]=p;
}
l=k*n+k;
a[l]=1.0/a[l];
for（j=0;j<=n-1;j++）
if（j!
=k）
{u=k*n+j;a[u]=a[u]*a[l];}
for（i=0;i<=n-1;i++）
if（i!
=k）
for（j=0;j<=n-1;j++）
if（j!
=k）
{u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
}
for（i=0;i<=n-1;i++）
if（i!
=k）
{u=i*n+k;a[u]=-a[u]*a[l];}
}
for（k=n-1;k>=0;k--）
{if（js[k]!
=k）
for（j=0;j<=n-1;j++）
{u=k*n+j;v=js[k]*n+j;
p=a[u];a[u]=a[v];a[v]=p;
}
if（is[k]!
=k）
for（i=0;i<=n-1;i++）
{u=i*n+k;v=i*n+is[k];
p=a[u];a[u]=a[v];a[v]=p;
}
}
for（i=0;i<=n-1;i++）
for（j=0;j<=n-1;j++）
*（A+i*n1+j）=*（a+i*n+j）;
free（is）;free（js）;
free（a）;
return
（1）;
}
//曲面拟合
//m为正交多项式系最高次值，array为拟合系数矩阵
//（u,v）为待拟和点
doubleqmlh（intm,doubleu,doublev,doublearray[]）
{
intr,s;
doublesum=0.0;
for（s=0;s{
for（r=0;r{
sum+=array[r*m+s]*pow（u,r）*pow（v,s）;
}
}
returnsum;
}
//Gauss列主元素消去法解线性方程组
//a,系数矩阵;b，值向量，返回时存放解向量
intagaus（a,b,n）
intn;
doublea[],b[];
{int*js,l,k,i,j,is,p,q;
doubled,t;
js=malloc（n*sizeof（int））;
l=1;
for（k=0;k<=n-2;k++）
{d=0.0;
for（i=k;i<=n-1;i++）
for（j=k;j<=n-1;j++）
{t=fabs（a[i*n+j]）;
if（t>d）{d=t;js[k]=j;is=i;}
}
if（d+1.0==1.0）l=0;
else
{if（js[k]!
=k）
for（i=0;i<=n-1;i++）
{p=i*n+k;q=i*n+js[k];
t=a[p];a[p]=a[q];a[q]=t;
}
if（is!
=k）
{for（j=k;j<=n-1;j++）
{p=k*n+j;q=is*n+j;
t=a[p];a[p]=a[q];a[q]=t;
}
t=b[k];b[k]=b[is];b[is]=t;
}
}
if（l==0）
{free（js）;printf（"fail\n"）;
return（0）;
}
d=a[k*n+k];
for（j=k+1;j<=n-1;j++）
{p=k*n+j;a[p]=a[p]/d;}
b[k]=b[k]/d;
for（i=k+1;i<=n-1;i++）
{for（j=k+1;j<=n-1;j++）
{p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[（n-1）*n+n-1];
if（fabs（d）+1.0==1.0）
{free（js）;printf（"fail\n"）;
return（0）;
}
b[n-1]=b[n-1]/d;
for（i=n-2;i>=0;i--）
{t=0.0;
for（j=i+1;j<=n-1;j++）
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for（k=n-1;k>=0;k--）
if（js[k]!
=k）
{t=b[k];b[k]=b[js[k]];b[js[k]]=t;}
free（js）;
return
（1）;
}
main（）
{
inti,j,r,s,b1,b2,k,ss,kk;
doublexx[11]={0.},yy[21]={0.},x[6]={0.},y[6]={0.},xin[8]={0.},yin[5]={0.};
doublea[4][4]={0.},b[4]={0.},h[4]={0.},H[4]={0.},g[4]={0.};
doubleC[11][11]={0.},B[11][11]={0.},BZ[11][11]={0.},BT[11][11]={0.};
doubleBTT[11][11]={0.},BU[11][21]={0.};
doubleG[21][11]={0.},GZ[11][21]={0.},GT[11][11]={0.};
doubleGTT[21][11]={0.},U[11][21]={0.};
doublebk,eps,u,v,uu,vv,w,delta,t1,t2,sigama,lh,llh;
doubleX[4]={0.5,1.5,-0.5,1.5};//初始迭代向量
doublez[6][6]={{-0.5,-0.42,-0.18,0.22,0.78,1.5},
{-0.34,-0.5,-0.5,-0.34,-0.02,0.46},{0.14,-0.26,-0.5,-0.58,-0.5,-0.26},
{0.94,0.3,-0.18,-0.5,-0.66,-0.66},{2.06,1.18,0.46,-0.1,-0.5,-0.74},{3.5,2.38,1.42,0.62,-0.02,-0.5}};//插值表
FILE*fp;
if（（fp=fopen（"SY0604229.txt","w+"））==NULL）//建立一个新文件存放输出结果
{
printf（"no\n"）;
}
bk=1.0;eps=1e-12;delta=10e-7;t1=0.;t2=0.;sigama=0.;lh=0.;
for（p1=0;p1<11;p1++）xx[p1]=0.08*p1;
for（q1=0;q1<21;q1++）yy[q1]=0.5+0.05*q1;
for（i=0;i<=5;i++）x[i]=0.4*i;
for（j=0;j<=5;j++）y[j]=0.2*j;
for（p1=0;p1<11;p1++）
{
for（q1=0;q1<21;q1++）
{
do
{
dnewta1（X,a）;
dnewta2（X,b,xx,yy）;
if（agaus（a,b,4）!
=0）
{
for（i=0;i<4;i++）
{
h[i]=X[i];
X[i]=X[i]+b[i];
H[i]=X[i];
g[i]=H[i]-h[i];
}
}
t1=cond（H,4）;
t2=cond（g,4）;
bk=t2/t1;//精度计算
}
while（bk>=eps）;//达到精度退出循环
u=X[1];v=X[0];//得出插值坐标
w=cz（x,y,z,u,v）;//对给定点进行插值计算
U[p1][q1]=w;
fprintf（fp,"x=%13.12ey=%13.12ef（x,y）=%13.12e\n",xx[p1],yy[q1],w）;
//将插值结果输出到文件SY0604229.txt
}
printf（"\n\n"）;
}
fprintf（fp,"\n\n\n"）;
do
{
for（k=1;k<11;k++）
{
for（j=0;j<=k;j++）
{
for（i=0;i<11;i++）
{
x[i]=0.08*i;//待拟合点的x坐标向量
B[i][j]=pow（x[i],j）;//x的各阶正交函数的幂次矩阵
}
}
for（b2=0;b2<=k;b2++）
{
for（b1=0;b1<21;b1++）
{
y[b1]=0.5+0.05*b1;//待拟合点的y坐标向量
G[b1][b2]=pow（y[b1],b2）;//y的各阶正交函数的幂次矩阵
}
}
zhuanzhi（B,BZ,11,11）;//对B矩阵进行转置，结果存放在BZ矩阵中
brmul（BZ,B,11,11,11,BT）;//将转置矩阵和B矩阵相乘，结果存放在BT矩阵中
ss=brinv（BT,k+1,11）;//对BT矩阵进行求逆，结果将原矩阵覆盖
if（ss!
=0）
{
brmul（BT,BZ,11,11,11,BTT）;
//将BT矩阵和BZ矩阵相乘，结果存放在BTT矩阵中
}
zhuanzhi（G,GZ,21,11）;//对G矩阵进行转置，结果存放在GZ矩阵中
brmul（GZ,G,11,21,11,GT）;//将转置矩阵和G矩阵相乘，结果存放在GT矩阵中
kk=brinv（GT,k+1,11）;//对GT矩阵进行求逆，结果将原矩阵覆盖
if（kk!
=0）
{
brmul（G,GT,21,11,11,GTT）;
//将G矩阵和GT矩阵相乘，结果存放在GTT矩阵中
}
brmul（BTT,U,11,11,21,BU）;
//将BTT矩阵和U矩阵相乘，结果存放在BU矩阵中
brmul（BU,GTT,11,21,11,C）;
//将BU矩阵和GTT矩阵相乘，结果存放在C矩阵中
for（i=0;i<11;i++）
{
for（j=0;j<21;j++）
{
u=0.08*i;
v=0.5+0.05*j;
lh=qmlh（11,u,v,C）;//对给定点进行曲面拟和
sigama+=（lh-U[i][j]）*（lh-U[i][j]）;//计算误差
}
}
}
}
while（sigama>=1e-7）;//达到精度退出循环
for（i=0;i<=k;i++）
{
for（j=0;j<=k;j++）
{
fprintf（fp,"%13.12e",C[i][j]）；
//将拟合系数矩阵输出到文件SY0604229.txt
}
fprintf（fp,"\n"）;
}
fprintf（fp,"LH=%13.12esigama=%13.12e\n",lh,sigama）;
//将误差输出到文件SY0604229.txt
for（i=0;i<6;i++）
{
for（j=0;j<6;j++）
{
fprintf（fp,"C[%d][%d]=%13.12e\n",i,j,C[i][j]）;
}
}
fprintf（fp,"\n"）;
bk=0.0;
for（p1=0;p1<8;p1++）
xin[p1]=0.1*p1+0.1;//初始化x*向量
for（q1=0;q1<5;q1++）
yin[q1]=0.7+0.2*q1;//初始化y*向量
for（i=0;i<=5;i++）
x[i]=0.4*i;//初始化插值节点的x坐标
for（j=0;j<=5;j++）
y[j]=0.2*j;//初始化插值节点的y坐标
for（p1=0;p1<8;p1++）
{
for（q1=0;q1<5;q1++）
{
do
{
dnewta1（X,a）;
dnewta2（X,b,xin,yin）;
if（agaus（a,b,4）!
=0）
{
for（i=0;i<4;i++）
{
h[i]=X[i];
X[i]=X[i]+b[i];
H[i]=X[i];
g[i]=H[i]-h[i];
}
}
t1=cond（H,4）;
t2=cond（g,4）;
bk=t2/t1;//精度计算
}
while（bk>=eps）;//达到精度退出循环
u=X[1];v=X[
升级会员 