北航数值分析大作业(三)

北航数值分析大作业（三）（看到那个收10个财富值的就不舒服，共享的精髓就在于免费）一、题目关于x，y，t，u，v，w的下列方程组0.5cost+u+v+w-x=2.67t+0.5sinu+v+w-y=1.070.5t+u+cosv+w-x=3.74t+0.5u+v+sinw-y=0.79以及关于z，t，u的下列二维数表zut00.40.81.21.620-0.5-0.340.140.942.063.50.2-0.42-0.5-0.260.31.182.380.4-0.18-0.5-0.5-0.180.461.420.60.22-0.34-0.58-0.5-0.10.620.80.78-0.02-0.5-0.66-0.5-0.021.01.50.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}上的一个近似表达式,0(,)krsrsrspxycxy要求p(x,y)一最小的k值达到以下的精度10202700((,)(,))10ijijijfxypxy其中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.将0.08(0,1,,10)ixii和0.50.05(0,1,,20)jyjj代入非线性方程组中，用牛顿法解出it和ju；2.以采取分片二次插值，选择（m，n）满足,2322mimhhtttm,2322njnuuun如果12ihtt或42ihtt，则m=1或4；如果12juu或42juu，则n=1或n=4。选择(,)(1,,1;1,,1)krtukmmmrnnn为插值节点，相应的Lagrange形式的插值多项式为112211(,)()()(,)mnkrkrkmrnptultluftu其中11()mwkwmkwwkttlttt(k=m-1,m,m+1)11()nwrwnrwwryyluyy(r=n-1,n,n+1)并将it和ju代入22(,)ptu，便得到了数表,,(,)ijijxyfxy。3.进行曲面拟和系数矩阵[]rscC，11()()TTTCBBBUGGG其中0011101011[()]1kkrikxxxxxxxB，0011101011[()]1kksjkyyyyGyyy[(,)]ijfxyUk从0逐渐增大，直到710，便得到了要求精度的系数rsc。三、全部源程序//在visualStudio2010的c语言环境下编译通过#includemath.h#includestdio.h#includeconio.h#definen4staticdoublex[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;in;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;in;i++){x_[i]=x[i]=1;}for(i=0;i2000;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;jn;j++){x[j]+=x_[j];}}cout迭代2000次没有带到精度要求endl;return;}voidSolve_tu()//求解对应的t、u子程序{inti,j;doublex_x[11],y_y[21];for(i=0;i11;i++)x_x[i]=0.08*i;for(j=0;j21;j++)y_y[j]=0.5+0.05*j;for(i=0;i11;i++)for(j=0;j21;j++){NEWTON(x_x[i],y_y[j]);t[i][j]=x[0];u[i][j]=x[1];}for(i=0;i8;i++)x_x[i]=0.1*(i+1);for(j=0;j5;j++)y_y[j]=0.5+0.2*(j+1);for(i=0;i8;i++)for(j=0;j5;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;ila;i++){for(j=0;jlb;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;k3;k++){temp1=1;temp2=1;for(m=0;m3;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;k3;k++)for(m=0;m3;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