北航数值分析大题目三.docx
- 文档编号:24475998
- 上传时间:2023-05-27
- 格式:DOCX
- 页数:30
- 大小:73.90KB
北航数值分析大题目三.docx
《北航数值分析大题目三.docx》由会员分享,可在线阅读,更多相关《北航数值分析大题目三.docx(30页珍藏版)》请在冰豆网上搜索。
北航数值分析大题目三
北航数值分析大作业(三)
一、题目
关于x,y,t,u,v,w的下列方程组
0.5cost+u+v+w-x=2.67
t+0.5sinu+v+w-y=1.07
0.5t+u+cosv+w-x=3.74
t+0.5u+v+sinw-y=0.79
以及关于z,t,u的下列二维数表
zu
t
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
确定了一个二元函数z=f(x,y)。
1.试用数值方法求出f(x,y)在区域D={(x,y)︱0≤x≤0.8,0.5≤y≤1.5}上的一个近似表达式
要求p(x,y)一最小的k值达到以下的精度
其中xi=0.08i,yj=0.5+0.05j。
2.计算f(xi*,yj*),p(xi*,yj*)(i=1,2,…,8;j=1,2,…,5)的值,以观察p(x,y)逼近f(x,y)的效果,其中xi*=0.1i,yj*=0.5+0.2j。
二、算法设计方案
1.将
和
代入非线性方程组中,用牛顿法解出
和
;
2.以采取分片二次插值,选择(m,n)满足
如果
或
,则m=1或4;如果
或
,则n=1或n=4。
选择
为插值节点,相应的Lagrange形式的插值多项式为
其中
(k=m-1,m,m+1)
(r=n-1,n,n+1)
并将
和
代入
,便得到了数表
。
3.进行曲面拟和系数矩阵
,
其中
,
k从0逐渐增大,直到
,便得到了要求精度的系数
。
三、全部源程序
//在visualStudio2010的c语言环境下编译通过
#include
#include
#include
#definen4
staticdoublex[n],x_[n],U[11][21],t[11][21],u[11][21],tt[8][5],uu[8][5],UU[8][5],c[9][9];
intk;
doubleMAX(doublea[n])//求数组中的最大值
{
inti;
doublemax;
max=fabs(a[0]);
for(i=0;i if(fabs(a[i])>max) max=fabs(a[i]); return(max); } voidDooLittle(doublea[n][n],doubleb[n])//选主元的DooLittle分解法求线性方程组 { inti,j,k,t,i_k,m[n]; doubleu[n][n],l[n][n],s[n],y[n]; doubleu_x,l_u,max,temp; for(k=1;k<=n;k++) { for(i=k;i<=n;i++) { l_u=0; for(t=1;t<=k-1;t++) l_u=l_u+l[i-1][t-1]*u[t-1][k-1]; s[i-1]=a[i-1][k-1]-l_u; } max=fabs(s[k-1]); i_k=k; m[k-1]=k; for(i=k;i<=n;i++) { if(fabs(s[i-1])>max) { max=fabs(s[i-1]); i_k=i; m[k-1]=i_k; } } if(i_k! =k) { for(t=1;t<=k-1;t++) { temp=l[k-1][t-1]; l[k-1][t-1]=l[i_k-1][t-1]; l[i_k-1][t-1]=temp; } for(t=k;t<=n;t++) { temp=a[k-1][t-1]; a[k-1][t-1]=a[i_k-1][t-1]; a[i_k-1][t-1]=temp; } temp=s[k-1]; s[k-1]=s[i_k-1]; s[i_k-1]=temp; } u[k-1][k-1]=s[k-1]; for(j=k+1;j<=n;j++) { l_u=0; for(t=1;t<=k-1;t++) l_u=l_u+l[k-1][t-1]*u[t-1][j-1]; u[k-1][j-1]=a[k-1][j-1]-l_u; } for(i=k+1;i<=n;i++) l[i-1][k-1]=s[i-1]/u[k-1][k-1]; } for(k=1;k<=n-1;k++) { temp=b[k-1]; b[k-1]=b[m[k-1]-1]; b[m[k-1]-1]=temp; } y[0]=b[0]; for(i=2;i<=n;i++) { l_u=0; for(t=1;t<=i-1;t++) l_u=l_u+l[i-1][t-1]*y[t-1]; y[i-1]=b[i-1]-l_u; } x_[n-1]=y[n-1]/u[n-1][n-1]; for(i=n-1;i>=1;i--) { u_x=0; for(t=i+1;t<=n;t++) u_x=u_x+u[i-1][t-1]*x_[t-1]; x_[i-1]=(y[i-1]-u_x)/u[i-1][i-1]; } } voidNEWTON(doublex_x,doubley_y)//牛顿法求解非线性方程组子程序 { doubleF[n],F_1[n][n]; inti,j; for(i=0;i { x_[i]=x[i]=1; } for(i=0;i<2000;i++) { F[0]=(2.67-0.5*cos(x[0])-x[1]-x[2]-x[3]+x_x); F[1]=(1.07-x[0]-0.5*sin(x[1])-x[2]-x[3]+y_y); F[2]=(3.74-0.5*x[0]-x[1]-cos(x[2])-x[3]+x_x); F[3]=(0.79-x[0]-0.5*x[1]-x[2]-sin(x[3])+y_y); F_1[0][0]=-0.5*sin(x[0]); F_1[0][1]=1 ;F_1[0][2]=1; F_1[0][3]=1; F_1[1][0]=1; F_1[1][1]=0.5*cos(x[1]); F_1[1][2]=1; F_1[1][3]=1; F_1[2][0]=0.5; F_1[2][1]=1; F_1[2][2]=-sin(x[2]);F_1[2][3]=1; F_1[3][0]=1; F_1[3][1]=0.5; F_1[3][2]=1; F_1[3][3]=cos(x[3]); DooLittle(F_1,F); if(MAX(x_)/MAX(x)<=1.0e-12) return; for(j=0;j { x[j]+=x_[j]; } } cout<<"迭代2000次没有带到精度要求"< return; } voidSolve_tu()//求解对应的t、u子程序 { inti,j; doublex_x[11],y_y[21]; for(i=0;i<11;i++) x_x[i]=0.08*i; for(j=0;j<21;j++) y_y[j]=0.5+0.05*j; for(i=0;i<11;i++) for(j=0;j<21;j++) { NEWTON(x_x[i],y_y[j]); t[i][j]=x[0]; u[i][j]=x[1]; } for(i=0;i<8;i++) x_x[i]=0.1*(i+1); for(j=0;j<5;j++) y_y[j]=0.5+0.2*(j+1); for(i=0;i<8;i++) for(j=0;j<5;j++) { NEWTON(x_x[i],y_y[j]); tt[i][j]=x[0]; uu[i][j]=x[1]; } } voidLagrange(double*t,double*u,double*U,intla,intlb,intflag)//分片二次代数插值 { inti,j,k,m,p,c,d,q; doublea[11][21],b[11][21],Z[6][6],temp1,temp2,L1[3],L2[3]; Z[0][0]=-0.5;Z[0][1]=-0.34;Z[0][2]=0.14;Z[0][3]=0.94;Z[0][4]=2.06;Z[0][5]=3.5; Z[1][0]=-0.42;Z[1][1]=-0.5;Z[1][2]=-0.26;Z[1][3]=0.3;Z[1][4]=1.18;Z[1][5]=2.38; Z[2][0]=-0.18;Z[2][1]=-0.5;Z[2][2]=-0.5;Z[2][3]=-0.18;Z[2][4]=0.46;Z[2][5]=1.42; Z[3][0]=0.22;Z[3][1]=-0.34;Z[3][2]=-0.58;Z[3][3]=-0.5;Z[3][4]=-0.1;Z[3][5]=0.62; Z[4][0]=0.78;Z[4][1]=-0.02;Z[4][2]=-0.5;Z[4][3]=-0.66;Z[4][4]=-0.5;Z[4][5]=-0.02; Z[5][0]=1.5;Z[5][1]=0.46;Z[5][2]=-0.26;Z[5][3]=-0.66;Z[5][4]=-0.74;Z[5][5]=-0.5; for(i=0;i { for(j=0;j { c=int(*(u+i*lb+j)/0.4); d=int((*(u+i*lb+j)-c*0.4)/(0.5*0.4)); a[i][j]=(c+d)*0.4; c=int(*(t+i*lb+j)/0.2); d=int((*(t+i*lb+j)-c*0.2)/(0.5*0.2)); b[i][j]=(c+d)*0.2; if(a[i][j]<0.4)a[i][j]=0.4; if(a[i][j]>1.6)a[i][j]=1.6; if(b[i][j]<0.2)b[i][j]=0.2; if(b[i][j]>0.8)b[i][j]=0.8; for(k=0;k<3;k++) { temp1=1;temp2=1; for(m=0;m<3;m++) if(m! =k) { temp1*=(*(t+i*lb+j)-(b[i][j]+(m-1)*0.2))/(b[i][j]+(k-1)*0.2-(b[i][j]+(m-1)*0.2)); temp2*=(*(u+i*lb+j)-(a[i][j]+(m-1)*0.4))/(a[i][j]+(k-1)*0.4-(a[i][j]+(m-1)*0.4)); } L1[k]=temp1; L2[k]=temp2; } temp1=0; for(k=0;k<3;k++) for(m=0;m<3;m++) { p=int(b[i][j]/0.2)-1+k; q=int(a[i][j]/0.4)-1+m; temp1+=L1[k]*L2[m]*Z[p][q]; } *(U+i*lb+j)=temp1; if(flag) { printf("%s%2d%s%.2f%s%2d%s%.2f%s%19.11e%s","x(",i,")=",0.08*i,",y(",j,")=",0.5+0.05*j,",",*(U+i*lb+j),";"); if(! fmod(j+1,2)) cout< if(j==lb-1) cout< } } } } voidMulti(double*a,double*b,double*c,intla,intlb,intlc,intr,ints,intt) //求r*s阶矩阵A与s*t阶矩阵B的乘积矩阵C { inti,j,k; for(i=0;i for(j=0;j { *(c+i*lc+j)=0; for(k=0;k } } doubleInverse(double*a,double*b,intla,intlb,intN)//求n阶方阵A的逆矩阵B { inti,j,k; doubletemp; for(i=0;i for(j=0;j if(i==j) *(b+i*lb+j)=1; else *(b+i*lb+j)=0; for(k=0;k j=k; for(i=k+1;i if(fabs(*(a+i*la+k))>fabs(*(a+j*la+k)))j=i; if(j! =k) for(i=0;i { temp=*(a+j*la+i); *(a+j*la+i)=*(a+k*la+i); *(a+k*la+i)=temp; temp=*(b+j*lb+i); *(b+j*lb+i)=*(b+k*lb+i); *(b+k*lb+i)=temp; } if(*(a+k*la+k)==0) return0; if((temp=*(a+k*la+k))! =1) for(i=0;i *(a+k*la+i)/=temp; *(b+k*lb+i)/=temp; } for(i=0;i if((*(a+i*la+k)! =0)&&(i! =k)){ temp=*(a+i*la+k); for(j=0;j *(a+i*la+j)-=temp*(*(a+k*la+j)); *(b+i*lb+j)-=temp*(*(b+k*lb+j)); } } } return0; } voidsolve_C() { inti,j,r,s; doublet1[21][21],t2[21][21],t3[21][21],d[9][9],e[9][9]; doubleB[11][9],B_T[9][11],G[21][9],G_T[9][21],P[11][21]; doubletemp,FangCha; for(i=0;i<9;i++) { for(j=0;j<11;j++) { B[j][i]=pow(0.08*j,i); B_T[i][j]=pow(0.08*j,i); } for(j=0;j<21;j++) { G[j][i]=pow(0.5+0.05*j,i); G_T[i][j]=pow(0.5+0.05*j,i); } } for(k=0;k<9;k++) { FangCha=0; Multi(B_T[0],B[0],t1[0],11,9,21,k+1,11,k+1); Inverse(t1[0],c[0],21,9,k+1); Multi(e[0],c[0],d[0],9,9,9,k+1,k+1,k+1); Multi(c[0],B_T[0],t1[0],9,11,21,k+1,k+1,11); Multi(t1[0],U[0],t2[0],21,21,21,k+1,11,21); Multi(G_T[0],G[0],t1[0],21,9,21,k+1,21,k+1); Inverse(t1[0],c[0],21,9,k+1); Multi(G[0],c[0],t3[0],9,9,21,21,k+1,k+1); Multi(t2[0],t3[0],c[0],21,21,9,k+1,21,k+1); for(i=0;i<11;i++) for(j=0;j<21;j++) { temp=0; for(r=0;r for(s=0;s temp+=c[r][s]*B[i][r]*G[j][s]; P[i][j]=temp; FangCha+=(U[i][j]-temp)*(U[i][j]-temp); } printf("%s%d%s%.11e%s","k=",k,";Sigma=",FangCa,";\n"); if(FangCha<=1.0e-7) { printf("%s%d%s%.11e%s","达到精度要求时: k=",k,";sigma=",FangCha,";\n系数c(r,s)如下: \n"); for(i=0;i { for(j=0;j { printf("%s%d%s%d%s%19.11e%s","c(",i,",",j,")=",c[i][j],";"); if(j==(k-1)/2) cout< } cout< } cout< return; } } cout<<"经过8次拟合没有达到所需精度;"< return; } voidCheck()//验证子程序 { inti,j,r,s; doubleB[8][9],G[5][9],temp; Lagrange(tt[0],uu[0],UU[0],8,5,0); for(i=0;i { for(j=0;j<8;j++) { B[j][i]=pow(0.1*(j+1),i); } for(j=0;j<5;j++) { G[j][i]=pow(0.5+0.2*(j+1),i); } } for(i=0;i<8;i++) for(j=0;j<5;j++) { temp=0; for(r=0;r for(s=0;s temp+=c[r][s]*B[i][r]*G[j][s]; printf("%s%d%s%.2f%s%d%s%.2f","x*(",i+1,")=",0.1*(i+1),",y*(",j+1,")=",0.5+0.2*(j+1)); printf("%s%19.11e%s%19.11e%s",",f(x*,y*)=",UU[i][j],",p(x*,y*)=",temp,";\n"); } } voidmain() { voidNEWTON(double,double); voidDooLittle(double**,double*); voidSolve_tu(); voidsolve_C(); voidCheck(); voidLagrange(double*t,double*u,double*U,intla,intlb,int); Solve_tu(); Lagrange(t[0],u[0],U[0],11,21,1); solve_C(); Check(); getchar(); } 四、程序运行结果: 1.数表: xi,yj,f(xi,yj) x(0)=0.00,y(0)=0.50,4.46504018481e-001;x(0)=0.00,y (1)=0.55,3.24683262928e-001; x(0)=0.00,y (2)=0.60,2.10159686683e-001;x(0)=0.00,y(3)=0.65,1.03043608316e-001; x(0)=0.00,y(4)=0.70,3.40189556268e-003;x(0)=0.00,y(5)=0.75,-8.87358136380e-002; x(0)=0.00,y(6)=0.80,-1.73371632750e-001;x(0)=0.00,y(7)=0.85,-2.50534611467e-001; x(0)=0.00,y(8)=0.90,-3.20276506388e-001;x(0)=0.00,y(9)=0.95,-3.82668066110e-001; x(0)=0.00,y(10)=1.00,-4.37795766738e-001;x(0)=0.00,y(11)=1.05,-4.85758941444e-001; x(0)=0.00,y(12)=1.10,-5.26667254884e-001;x(0)=0.00,y(13)=1.15,-5.60638479797e-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 题目