北航数值分析大作业3文档格式.docx
- 文档编号:16681263
- 上传时间:2022-11-25
- 格式:DOCX
- 页数:30
- 大小:111.52KB
北航数值分析大作业3文档格式.docx
《北航数值分析大作业3文档格式.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业3文档格式.docx(30页珍藏版)》请在冰豆网上搜索。
,若
,则选
;
若
,等上述方法一样。
其中:
4.二元函数拟合
已通过求解非线性方程组和插值得到的11×
21的((x,y),z)二元矩形数表为基准,进行给定不同次数的二元函数最小二乘法拟合。
拟合过程选取正交函数系作为基函数,即拟合函数可表示为:
,其中
和
分别为次数为r和s的正交多项式系。
将其写为矩阵形式为:
同理:
所以:
其中
利用x和y的正交函数系同其幂函数系数阵以及
(2)中所求系数矩阵,最终得到所求系数矩阵,即:
该子程序中的矩阵乘积运算调用了通用的矩阵相乘的子程序。
二、全部源程序
#include"
stdio.h"
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;
for(j=2;
j<
j++)
if(v>
(y[j]-f/2)&
v<
=(y[j]+f/2))
k2=j;
for(ii=k1-1;
ii<
k1+2;
ii++)
bb[ii]=0.0;
for(jj=k2-1;
jj<
k2+2;
jj++)
{
k=ii*6+jj;
hh=z[k];
for(kk=k2-1;
kk<
kk++)
if(kk!
=jj)
hh*=(v-y[kk])/(y[jj]-y[kk]);
bb[ii]=bb[ii]+hh;
}
w=0.0;
=k1+1;
hh=bb[ii];
for(jj=k1-1;
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;
//求向量的无穷范数
//array为向量地址,n为向量维数
doublecond(doublearray[],intn)
inti;
doublefan=1e-32;
for(i=0;
n;
if(fabs(array[i])>
fan)
fan=fabs(array[i]);
returnfan;
}//矩阵转置
//a为待转置矩阵,维数为m*n,转置后存放在矩阵b中
zhuanzhi(doublea[],doubleb[],intm,intn)
inti,j;
m;
for(j=0;
b[m*j+i]=a[i*n+j];
}
}//矩阵相乘
//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];
//对矩阵部分求逆
//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));
*(a+i*n+j)=*(A+i*n1+j);
for(k=0;
k<
k++)
{d=0.0;
for(i=k;
for(j=k;
{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)
{u=k*n+j;
v=is[k]*n+j;
p=a[u];
a[u]=a[v];
a[v]=p;
if(js[k]!
{u=i*n+k;
v=i*n+js[k];
l=k*n+k;
a[l]=1.0/a[l];
if(j!
a[u]=a[u]*a[l];
if(i!
{u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
a[u]=-a[u]*a[l];
for(k=n-1;
k>
=0;
k--)
{if(js[k]!
v=js[k]*n+j;
v=i*n+is[k];
*(A+i*n1+j)=*(a+i*n+j);
free(is);
free(a);
return
(1);
//曲面拟合
//m为正交多项式系最高次值,array为拟合系数矩阵
//(u,v)为待拟和点
doubleqmlh(intm,doubleu,doublev,doublearray[])
intr,s;
doublesum=0.0;
for(s=0;
s<
s++)
for(r=0;
r<
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;
l=1;
k<
=n-2;
k++)
{t=fabs(a[i*n+j]);
if(t>
d){d=t;
is=i;
if(d+1.0==1.0)l=0;
{p=i*n+k;
q=i*n+js[k];
t=a[p];
a[p]=a[q];
a[q]=t;
if(is!
{for(j=k;
{p=k*n+j;
q=is*n+j;
t=b[k];
b[k]=b[is];
b[is]=t;
if(l==0)
{free(js);
printf("
fail\n"
d=a[k*n+k];
for(j=k+1;
a[p]=a[p]/d;
b[k]=b[k]/d;
for(i=k+1;
{for(j=k+1;
{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)
b[n-1]=b[n-1]/d;
for(i=n-2;
i>
i--)
{t=0.0;
for(j=i+1;
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
js[n-1]=n-1;
k>
k--)
{t=b[k];
b[k]=b[js[k]];
b[js[k]]=t;
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)//建立一个新文件存放输出结果
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;
=5;
i++)x[i]=0.4*i;
for(j=0;
j++)y[j]=0.2*j;
p1++)
for(q1=0;
q1++)
do
{
dnewta1(X,a);
dnewta2(X,b,xx,yy);
if(agaus(a,b,4)!
=0)
{
for(i=0;
4;
{
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
\n\n"
fprintf(fp,"
\n\n\n"
do
for(k=1;
for(j=0;
=k;
for(i=0;
x[i]=0.08*i;
//待拟合点的x坐标向量
B[i][j]=pow(x[i],j);
//x的各阶正交函数的幂次矩阵
for(b2=0;
b2<
b2++)
for(b1=0;
b1<
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!
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!
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;
for(j=0;
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);
%13.12e"
C[i][j]);
//将拟合系数矩阵输出到文件SY0604229.txt
\n"
LH=%13.12esigama=%13.12e\n"
lh,sigama);
//将误差输出到文件SY0604229.txt
6;
C[%d][%d]=%13.12e\n"
i,j,C[i][j]);
bk=0.0;
8;
xin[p1]=0.1*p1+0.1;
//初始化x*向量
5;
q1++)
yin[q1]=0.7+0.2*q1;
//初始化y*向量
i++)
x[i]=0.4*i;
//初始化插值节点的x坐标
j++)
y[j]=0.2*j;
//初始化插值节点的y坐标
dnewta2(X,b,xin,yin);
v=X[
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业