数值分析实习作业三

数值分析大作业第三题学号：姓名：学院：授课教师：吕淑娟2016.11.21题目：关于x,y,t,u,v,w的下列方程组（A.3）0.5cos2.670.5sin1.070.5cos3.740.5sin0.79tuvwxtuvwytuvwxtuvwy(A.3)以及关于z,t,u的下列二维数表zut00.40.81.21.62.00-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)在区域{(,)|00.8,0.51.5}Dxyxy上的一个近似表达式,0(,)krsrsrspxycxy要求(,)pxy以最小的k值达到以下的精度10202700((,)(,))10ijijijfxypxy其中，0.08,0.50.05ijxiyj。2、计算****(,),(,)ijijfxypxy（i=1,2,…，8；j=1,2,…，5）的值，以观察(,)pxy逼近(,)fxy的效果，其中，*ix=0.1i,*jy=0.5+0.2j。说明：1、用迭代方法求解非线性方程组时，要求近似解向量()kx满足()(1)()12||||/||||10kkkxxx2、作二元插值时，要使用分片二次代数插值。3、要由程序自动确定最小的k值。4、打印以下内容：（1）算法的设计方案；（2）全部源程序；（3）数表：,,ijxy(,)ijfxy（i=0,1,2,…,10；j=0,1,2,…,20）；（4）选择过程的,k值；（5）达到精度要求时的,k值以及(,)pxy中的系数rsc（r=0,1,…,k；s=0,1,…,k）；（6）数表：**,,ijxy****(,),(,)ijijfxypxy（i=1,2,…,8；j=1,2,…,5）；（7）讨论：计算过程中遇到的现象。5、采用f型输出,,ijxy**,ijxy的准确值，其余实型数采用e型输出并且至少显示12位有效数字。一、算法设计方案1、方程组(A.3)由四个方程组成，又为六元方程组，则给(𝑥,𝑦)赋一组值，则能解出一组（𝑡,𝑢,𝑣,𝑤）。依次代入𝑥𝑖=0.08𝑖，𝑦𝑗=0.5+0.05𝑗(其中𝑖=0,1,…,10;𝑗=0,1,…,20)，利用NEWTON迭代法解得对应的(𝑡𝑖𝑗,𝑢𝑖𝑗)。NEWTON迭代法具体算法见教材p90~91。2、利用二维数表进行分片二次代数插值，得到函数𝑧=𝑧(𝑡,𝑢)。对于区域D内任意(𝑥,𝑦)，代入方程组（A.3）解得(𝑡,𝑢)，再代入函数𝑧=𝑧(𝑡,𝑢)，此从(𝑥,𝑦)到𝑧的映射即为二元函数𝑧=𝑓(𝑥,𝑦)。分片二次代数插值具体算法见教材p101。3、将1中得到的(𝑡𝑖𝑗,𝑢𝑖𝑗)代入函数𝑧=𝑧(𝑡,𝑢)，得到𝑧𝑖𝑗=𝑓(𝑥𝑖,𝑦𝑗)。利用最小二乘法进行拟合，得到,0(,)krsrsrspxycxy，其中使𝑘依次取1,2,…直到使拟合函数满足10202700((,)(,))10ijijijfxypxy。则所取𝑘为满足条件的最小值。最小二乘法拟合具体算法见教材p142。4、同1、2步，代入*ix=0.1i,*jy=0.5+0.2j，（i=1,2,…，8；j=1,2,…，5），得到**(,)ijfxy；将**(,)ijxy代入(,)pxy，得到**(,)ijpxy。比较两函数值，观察(,)pxy逼近(,)fxy的效果。二、结果1、数表：,,ijxy(,)ijfxy（i=0,1,2,…,10；j=0,1,2,…,20）2、选择过程的,k值3、最终k应为5，此时crs矩阵如下4、数表：**,,ijxy****(,),(,)ijijfxypxy（i=1,2,…,8；j=1,2,…,5）三、遇到的问题从插值得到的f(x,y)函数开始（包括其在内），所得到计算结果与课本答案存在1e-9左右的差异。经验证，并不是迭代初值的选择所导致的。四、源程序#includestdafx.h#includestdio.h#includemath.h#includestdlib.h#defineN4#defineL11#defineM21int_tmain(intargc,_TCHAR*argv[]){voidntsl(doublex,doubley,doublea[N]);//Newton迭代解方程组函数voidgsslm(doublea[M][M],doubleb[M],doublex[M],intn);voidgssln(doublea[N][N],doubleb[N],doublex[N]);//列主元gauss迭代voidfcfz(doublea[N],doublef[N],doublex,doubley);//方程组赋值函数voidfcds(doublea[N],doubleff[N][N]);//方程导数矩阵函数doubleamul(doublea[N],doubleb[N]);//向量内积函数doublefpcz(doublet0[3],doubleu0[3],doublez0[3][3],doublet,doubleu);//二次分片插值函数voidatra(doublea[M][M],intn,intm,doubleat[M][M]);//矩阵转置voidaamu(doublea[M][M],intn,intm,ints,doubleb[M][M],doublec[M][M]);//矩阵相乘voidaqni(doublea[M][M],intn,doubleani[M][M]);//矩阵求逆doublepxy(doublex,doubley,doublec[M][M],intk);//求p(x,y)voidzxec(intl,intm,doubleu[L][M],doublex[L],doubley[M],doublec[M][M]);//最小二乘doubleztu(doublet,doubleu,doublez00[6][6]);//求z(t,u)doubletuvw[N]={1,1,1,1};//方程组迭代赋初值doublex;doubley;doublez[L][M];intk0;doublec[M][M]={0};doublex1[L];doubley1[M];doublez00[6][6]={{-0.5,-0.34,0.14,0.94,2.06,3.5},{-0.42,-0.5,-0.26,0.3,1.18,2.38},{-0.18,-0.5,-0.5,-0.18,0.46,1.42},{0.22,-0.34,-0.58,-0.5,-0.1,0.62},{0.78,-0.02,-0.5,-0.66,-0.5,-0.02},{1.5,0.46,-0.26,-0.66,-0.74,-0.5}};for(inti=0;iL;i++){for(intj=0;j7;j++){if(i==0){if(j==0)printf(x-y);printf(%1.2lf,0.5+0.05*j);if(j==6)printf(\n);}}for(intj=0;j7;j++){x=0.08*i;y=0.5+0.05*j;ntsl(x,y,tuvw);z[i][j]=ztu(tuvw[0],tuvw[1],z00);if(j==0)printf(%1.2lf,0.08*i);printf(%1.11e,z[i][j]);if(j==6)printf(\n);}if(i==10)printf(\n);}for(inti=0;iL;i++){for(intj=7;j14;j++){if(i==0){if(j==7)printf(x-y);printf(%1.2lf,0.5+0.05*j);if(j==13)printf(\n);}}for(intj=7;j14;j++){x=0.08*i;y=0.5+0.05*j;ntsl(x,y,tuvw);z[i][j]=ztu(tuvw[0],tuvw[1],z00);if(j==7)printf(%1.2lf,0.08*i);printf(%1.11e,z[i][j]);if(j==13)printf(\n);}if(i==10)printf(\n);}for(inti=0;iL;i++){for(intj=14;j21;j++){if(i==0){if(j==14)printf(x-y);printf(%1.2lf,0.5+0.05*j);if(j==20)printf(\n);}}for(intj=14;j21;j++){x=0.08*i;y=0.5+0.05*j;ntsl(x,y,tuvw);z[i][j]=ztu(tuvw[0],tuvw[1],z00);if(j==14)printf(%1.2lf,0.08*i);printf(%1.11e,z[i][j]);if(j==20)printf(\n);}if(i==10)printf(\n);}for(inti=0;iL;i++)x1[i]=0.08*i;for(inti=0;iM;i++)y1[i]=0.5+0.05*i;//最小二乘法拟合+循环确定满足条件k值for(intk=0;kL-1;k++){doublesum1[L][M]={0};doublesum2=0;zxec(k,k,z,x1,y1,c);for(inti=0;iL;i++){for(intj=0;jM;j++){for(intr=0;rk+1;r++){for(ints=0;sk+1;s++){if((r==0)&&(i==0))sum1[i][j]=sum1[i][j]+c[r][s]*1*pow(y1[j],s);elsesum1[i][j]=sum1[i][j]+c[r][s]*pow(x1[i],r)*pow(y1[j],s);}}sum2=sum2+pow((sum1[i][j]-z[i][j]),2);}}printf(k=%d,k);printf(σ=%1.11e\n,sum2);if(sum2(1e-7)){printf(k=%d\n,k);k0=k;break;}if(k==L-2)printf(拟合失败\n);}for(inti=0;i6;i++){for(intj=0;j6;j++){printf(%.11e,c[i][j]);if(j==5)printf(\n);}}for(inti=1;i9;i++)//求f(x*,y*){if(i==1)printf(x-y-f);for(intj=1;j6;j++){if(i==1){printf(%.2lf,0.5+0.2*j);if(j==5)printf(\n);}}printf(%.2lf,0.1*i);for(intj=1;j6;j++){x=0.1*i;y=0.5+0.2*j;ntsl(x,y,tuvw);printf(%.11e,ztu(tuvw[0],tuvw[1],z00));if(j==5)printf(\n);}}printf(\n);//求p(x*,y*)for(inti=1;i9;i++){if(i==1)printf(x-y-p);for(intj=1;j6;j++){if(i==1){printf(%.2lf,0.5+0.2*j);if(j==5)printf(\n);}}printf(%.2lf,0.