北航数值分析大作业3.docx
- 文档编号:3774984
- 上传时间:2022-11-25
- 格式:DOCX
- 页数:30
- 大小:111.52KB
北航数值分析大作业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;kk if(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[
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业
![提示](https://static.bdocx.com/images/bang_tan.gif)