北航数值分析2-二次插值与曲面拟合

1、数值分析2015/12/6北京航空航天大学数值分析大作业姓名：郝海东学号：ZY1506227准备工作Ø 算法设计分析一：题目要求出(x,y,z)三者之间的关系，但首先给出的测量点是(u,t,z)。所以首先为得到点集(x,y,z)，考虑在方程组A.3中带入每一对u,t，得到一个包含v,w,x,y四个变元的四元非线性方程组，此时采用非线性方程组的迭代法计算出(v,w,x,y)，从而得到相应(u,t,z)对应的(x,y,z)。分析二：题目未给出f(x,y)的精确表达式，却要求以一定的精度拟合原方程。鉴于此处的考量，决定以二元插值法得到的表达式作为f(x,y)的精确表达式，以x,y的k次多项

2、式（k未知，需要在程序中从1次到n次循环计算从而得到最佳的k值）拟合函数作为p(x,y)。分析三： 非线性方程组的求解有普通迭代法和牛顿迭代法等。在考虑这两种方法时，首先考虑到普通迭代收敛速度为一阶收敛，且要求 |F(X)|F < 1才能收敛。而在本例中， 1 1 -1 0 1 1 0 -1F(X) =-sinv 1 -1 0 |F(X)|F < 1不成立，所以普通迭代法无法收敛。所以应 1 cosw 0 -1该采用牛顿迭代法计算(v,w,x,y)。综上可得，解决本问题的步骤为首先用牛顿迭代法求解非线性方程组得到(u,t,z)对应的(x,y,z)点集，然后用分片二元函数插值法得到方

3、程的精确表达式f(x,y)，最后以xk,yk（k=1,2,）为基函数拟合得到p(x,y)。Ø 模块设计根据以上分析，设计以下模块：数据录入初始化模块void init()迭代计算(x,y,z)入口函数，迭代计算(v,w,x,y)的控制主体void iterCalXY()被调用函数，线性方程组求解得到v,wx,yvoid linearCal()二元函数插值二元函数二次插值表达式，返回f(x,y)double interpolation(double x, double y)返回朗格朗日系数给插值表达式double l(int k,double x,int flag)曲面拟合拟合函数表达

4、式，返回p(x,y)double surfaceFit(int k,double x, double y)矩阵相乘void trmul( double *a,double *b,int m,int n,int k,double *c)矩阵求逆int rinv(double *a,int n)计算k值int calK()结果输出void put(int k)Ø 数据结构设计本例中未涉及到复杂数据结构，(x,y,z)点集和非线性方程组F(X)采用二维矩阵存储。在曲面拟合时，各种矩阵的存储为一维数组，在进行行列的转换也较为方便。完整代码（编译环境windows10 + visual stu

5、dio2010）：/ algo2.cpp : 定义控制台应用程序的入口点。#include "stdafx.h"#include <math.h>#include <malloc.h>double iter45;double point77;double deltav,deltaw,deltax,deltay;double *c;#define print printf("-n");void init()iter00 = 1;iter01 = 1;iter02 = -1;iter03 = 0;iter10 = 1;iter11 =

6、 1;iter12 = 0;iter13 = -1;iter20 = 0;iter21 = 1;iter22 = -1;iter23 = 0;iter30 = 1;iter31 = 0;iter32 = 0;iter33 = -1;iter04 = 1;iter14 = 1;iter24 = 1;iter34 = 1;point00 = 0.0;point01 = 0;point02 = 0.4;point03 = 0.8;point04 = 1.2;point05 = 1.6;point06 = 2;point10 = 0.0;point11 = -0.5;point12 = -0.34;p

7、oint13 = 0.14;point14 = 0.94;point15 = 2.06;point16 = 3.5;point20 = 0.2;point21 = -0.42;point22 = -0.5;point23 = -0.26;point24 = 0.3;point25 = 1.18;point26 = 2.38;point30 = 0.4;point31 = -0.18;point32 = -0.5;point33 = -0.5;point34 = -0.18;point35 = 0.46;point36 = 1.42;point40 = 0.5;point41 = 0.22;po

8、int42 = -0.34;point43 = -0.58;point44 = -0.5;point45 = -0.1;point46 = 0.62;point50 = 0.8;point51 = 0.78;point52 = -0.02;point53 = -0.5;point54 = -0.66;point55 = -0.5;point56 = -0.02;point60 = 1.0;point61 = 1.5;point62 = 0.46;point63 = -0.26;point64 = -0.66;point65 = -0.74;point66 = -0.5;void linearC

9、al()double max;int k;double temp;for(int j = 0 ; j < 3 ; j+)max = 0;for(int i = j ; i < 4 ; i+)if(fabs(iterij) > max)max = fabs(iterij);k = i;for(int m = j ; m < 5 ; m+)temp = iterjm;iterjm = iterkm;iterkm = temp;for(int m = j+1; m < 4 ; m+)double kk = itermj / iterjj;for(int n = j ;

10、n < 5 ; n+)itermn = itermn - kk * iterjn;deltay = iter34 / iter33;deltax = (iter24 - deltay * iter23) / iter22;deltaw = (iter14 - deltay * iter13 - deltax * iter12) / iter11;deltav = (iter04 - deltay * iter03 - deltax * iter02 - deltaw * iter01) / iter00;/*for(int i = 0 ; i < 4 ; i+)printf(&qu

11、ot;%.3f %.3f %.3f %.3f %.3fn",iteri0,iteri1,iteri2,iteri3,iteri4);printf("n%.3f %.3f %.3f %.3fn",deltav,deltaw,deltax,deltay);*/void iterCalXY()double v=1,w=1,x=1,y=1;printf("非线性方程组迭代计算v,w,x,y：nn");for(int n = 1 ; n <= 6 ; n+)int k = 0;double t = pointn0;double u = point0

12、n;doiter00 = 1;iter01 = 1;iter02 = -1;iter03 = 0;iter10 = 1;iter11 = 1;iter12 = 0;iter13 = -1;iter20 = 0;iter21 = 1;iter22 = -1;iter23 = 0;iter30 = 1;iter31 = 0;iter32 = 0;iter33 = -1;iter04 = -(0.5 * cos(t) + u + v + w - x - 2.67);iter14 = -(t + 0.5 * sin(u) + v + w - y - 1.07);iter24 = -(0.5 * t +

13、 u + cos(v) + w - x - 3.74);iter34 = -(t + 0.5 * u + v + sin(w) - y - 0.79);iter20 = -sin(v);iter31 = cos(w);linearCal();v += deltav;w += deltaw;x += deltax;y += deltay;k+;while(fabs(deltav / v) > pow(10.0,-12) | fabs(deltaw / w) > pow(10.0,-12) | fabs(deltax / x) > pow(10.0,-12) | fabs(del

14、tay / y) > pow(10.0,-12);printf(" v=%.12e w=%.12e x=%.12e y=%.12e 迭代次数k=%d n",v,w,x,y,k);pointn0 = y;point0n = x;print;double l(int k,double x,int flag)double numerator=1;double denominator=1;if(flag = 0)for(int i = 0 ; i < 3 ; i+)if(k != i)numerator = numerator * (x - point0i+3);den

15、ominator = denominator * (point0k+3 - point0i+3);elsefor(int i = 0 ; i < 3 ; i+)if(k != i)numerator = numerator * (x - pointi+30);denominator = denominator * (pointk+30 - pointi+30);return numerator / denominator;double interpolation(double x, double y)double ret = 0;for(int k = 0 ; k < 3 ; k+

16、)for(int r = 0 ; r < 3 ; r+)ret = ret + l(k,x,0) * l(r,y,1) * pointk+3r+3;return ret;void trmul( double *a,double *b,int m,int n,int k,double *c)int i,j,l,u;for (i=0; i<=m-1; i+)for (j=0; j<=k-1; j+) u=i*k+j; cu=0.0;for (l=0; l<=n-1; l+)cu=cu+ai*n+l*bl*k+j;int rinv(double *a,int n)int *i

17、s,*js,i,j,k,l,u,v; double d,p; is=(int *)malloc(n*sizeof(int); js=(int *)malloc(n*sizeof(int); 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(al); if (p>d) d=p; isk=i; jsk=j; if (d+1.0=1.0) free(is); free(js); printf("err*not invn"); re

18、turn(0); if (isk!=k) for (j=0; j<=n-1; j+) u=k*n+j; v=isk*n+j; p=au; au=av; av=p; if (jsk!=k) for (i=0; i<=n-1; i+) u=i*n+k; v=i*n+jsk; p=au; au=av; av=p; l=k*n+k; al=1.0/al; for (j=0; j<=n-1; j+) if (j!=k) u=k*n+j; au=au*al; for (i=0; i<=n-1; i+) if (i!=k) for (j=0; j<=n-1; j+) if (j

19、!=k) u=i*n+j; au=au-ai*n+k*ak*n+j; for (i=0; i<=n-1; i+) if (i!=k) u=i*n+k; au=-au*al; for (k=n-1; k>=0; k-) if (jsk!=k) for (j=0; j<=n-1; j+) u=k*n+j; v=jsk*n+j; p=au; au=av; av=p; if (isk!=k) for (i=0; i<=n-1; i+) u=i*n+k; v=i*n+isk; p=au; au=av; av=p; free(is); free(js); return(1);dou

20、ble surfaceFit(int k,double x, double y)double *b= (double *)malloc(sizeof(double3) * (k+1);double *g= (double *)malloc(sizeof(double3) * (k+1);double *bt= (double *)malloc(sizeof(double3) * (k+1);double *gt= (double *)malloc(sizeof(double3) * (k+1);for(int j = 0 ; j <= k ; j+)for(int i = 0 ; i &

21、lt; 3 ; i+)bi*(k+1)+j = pow(point0i+3,j);gi*(k+1)+j = pow(pointi+30,j);btj*3+i = pow(point0i+3,j);gtj*3+i = pow(pointi+30,j);double *btb = (double *)malloc(sizeof(double) * (k+1) * (k+1);double *gtg = (double *)malloc(sizeof(double) * (k+1) * (k+1);trmul(bt,b,k+1,3,k+1,btb);trmul(gt,g,k+1,3,k+1,gtg)

22、;rinv(btb,k+1);rinv(gtg,k+1);double *u = (double *)malloc(sizeof(double) * 9);for(int i = 0 ; i < 3 ; i+)for(int j = 0 ; j < 3 ;j+)ui*3+j = pointi+3j+3;double *btbrbt = (double *)malloc(sizeof(double) * (k+1) * 3);trmul(btb,bt,k+1,k+1,3,btbrbt);double *btbrbtu = (double *)malloc(sizeof(double)

23、 * (k+1) * 3);trmul(btbrbt,u,k+1,3,3,btbrbtu);double *btbrbtug = (double *)malloc(sizeof(double) * (k+1) * (k+1);trmul(btbrbtu,g,k+1,3,(k+1),btbrbtug);c = (double *)malloc(sizeof(double) * (k+1) * (k+1);trmul(btbrbtug,gtg,k+1,k+1,k+1,c);double ret=0;for(int s = 0 ; s <= k ; s+ )for(int r = 0 ; r

24、<= k ;r+)ret = ret + cs*(k+1)+r * pow(x,s) * pow(y,r);return ret;int calK()int k = 1;double sigma = 0;putchar(10);print;printf("计算k值.n",k,sigma);dosigma = 0;for(int i = 0 ; i <= 10 ; i+)for(int j = 0 ; j <= 20 ; j+)sigma = sigma + pow(surfaceFit(k,0.08*i,0.5+0.05*j) - interpolatio

25、n(0.08*i,0.5+0.05*j),2);printf(" k=%d sigma=%.12en",k,sigma);k+;while(fabs(sigma) > pow(10.0,-7);k-;putchar(10);printf("当k=%d时,sigma=%.12e,满足要求。此时的c(r,s)如下所示:n",k,sigma);for(int r = 0 ; r <= k ; r+)for(int s = 0 ; s <= k ; s+)printf(" c(%d,%d)=%.12e ",r,s,cr*(k

26、+1)+s);putchar(10);print;putchar(10);putchar(10);return k;void put(int k)print;printf("当x(i)=0.08 * i,y(j)=0.5 + 0.05 * j 时，采用分片二次插值得到以下结果:nn");for(int i = 0 ; i <= 10 ; i+)for(int j = 0 ; j <= 20 ; j+)printf(" x(%d)=%f,y(%d)=%f,f(x,y)=%.12en",i,0.08*i,j,0.5+0.05*j,interpol

27、ation(0.08*i,0.5+0.05*j);print;putchar(10);print;printf("当x(i)=0.1 * i,y(j)=0.5 + 0.2 * j 时，采用分片二次插值和k次曲面拟合，分别得到以下结果:nn");for(int i = 0 ; i <= 8 ; i+)for(int j = 0 ; j <= 5 ; j+)printf(" x(%d)=%f,y(%d)=%f,f(x,y)=%.12e,p(x,y)=%.12en",i,0.1*i,j,0.5+0.2*j,interpolation(0.1*i,0

28、.5+0.2*j),surfaceFit(k,0.1*i,0.5+0.2*j);print;int _tmain(int argc, _TCHAR* argv)init();iterCalXY();int k = calK();put(k);return 0;结果输出非线性方程组迭代计算v,w,x,y： v=-1.402288635333e+000 w=1.218626459132e+000 x=-2.353662176202e+000 y=-1.253662176202e+000 迭代次数k=12 v=-6.866751502385e-001 w=1.226657014608e+000 x=

29、-1.239984846710e+000 y=-1.353089644759e-001 迭代次数k=9 v=-4.156943533737e-001 w=1.279071060422e+000 x=-5.460927959499e-001 y=5.520547524985e-001 迭代次数k=6 v=-3.049204401019e-001 w=1.399313225520e+000 x=6.318406636288e-002 y=9.904123284013e-001 迭代次数k=5 v=-1.852493640097e-002 w=1.580169272743e+000 x=8.3999

30、76910156e-001 y=1.791431137863e+000 迭代次数k=6 v=1.487987274582e-001 w=1.799347059348e+000 x=1.548296939740e+000 y=2.332794500219e+000 迭代次数k=5-计算k值. k=1 sigma=2.820715164735e+000 k=2 sigma=5.007125499241e-026当k=2时,sigma=5.007125499241e-026,满足要求。此时的c(r,s)如下所示: c(0,0)=-5.775556867930e-001 c(0,1)=-1.30627

31、7917977e-001 c(0,2)=2.366859028134e-001 c(1,0)=5.751212454893e-001 c(1,1)=-1.299373662970e+000 c(1,2)=3.035984844691e-001 c(2,0)=6.391262749418e-002 c(2,1)=2.186046543440e-001 c(2,2)=-5.107694856998e-002-当x(i)=0.08 * i,y(j)=0.5 + 0.05 * j 时，采用分片二次插值得到以下结果: x(0)=0.000000,y(0)=0.500000,f(x,y)=-5.83698

32、1069886e-001 x(0)=0.000000,y(1)=0.550000,f(x,y)=-5.778034866807e-001 x(0)=0.000000,y(2)=0.600000,f(x,y)=-5.707254368588e-001 x(0)=0.000000,y(3)=0.650000,f(x,y)=-5.624639575229e-001 x(0)=0.000000,y(4)=0.700000,f(x,y)=-5.530190486729e-001 x(0)=0.000000,y(5)=0.750000,f(x,y)=-5.423907103088e-001 x(0)=0.

33、000000,y(6)=0.800000,f(x,y)=-5.305789424306e-001 x(0)=0.000000,y(7)=0.850000,f(x,y)=-5.175837450384e-001 x(0)=0.000000,y(8)=0.900000,f(x,y)=-5.034051181321e-001 x(0)=0.000000,y(9)=0.950000,f(x,y)=-4.880430617118e-001 x(0)=0.000000,y(10)=1.000000,f(x,y)=-4.714975757774e-001 x(0)=0.000000,y(11)=1.0500

34、00,f(x,y)=-4.537686603289e-001 x(0)=0.000000,y(12)=1.100000,f(x,y)=-4.348563153663e-001 x(0)=0.000000,y(13)=1.150000,f(x,y)=-4.147605408897e-001 x(0)=0.000000,y(14)=1.200000,f(x,y)=-3.934813368990e-001 x(0)=0.000000,y(15)=1.250000,f(x,y)=-3.710187033943e-001 x(0)=0.000000,y(16)=1.300000,f(x,y)=-3.47

35、3726403755e-001 x(0)=0.000000,y(17)=1.350000,f(x,y)=-3.225431478426e-001 x(0)=0.000000,y(18)=1.400000,f(x,y)=-2.965302257956e-001 x(0)=0.000000,y(19)=1.450000,f(x,y)=-2.693338742346e-001 x(0)=0.000000,y(20)=1.500000,f(x,y)=-2.409540931595e-001 x(1)=0.080000,y(0)=0.500000,f(x,y)=-5.825645315867e-001

36、x(1)=0.080000,y(1)=0.550000,f(x,y)=-5.805395006613e-001 x(1)=0.080000,y(2)=0.600000,f(x,y)=-5.772112352904e-001 x(1)=0.080000,y(3)=0.650000,f(x,y)=-5.725797354741e-001 x(1)=0.080000,y(4)=0.700000,f(x,y)=-5.666450012122e-001 x(1)=0.080000,y(5)=0.750000,f(x,y)=-5.594070325048e-001 x(1)=0.080000,y(6)=0

37、.800000,f(x,y)=-5.508658293519e-001 x(1)=0.080000,y(7)=0.850000,f(x,y)=-5.410213917535e-001 x(1)=0.080000,y(8)=0.900000,f(x,y)=-5.298737197097e-001 x(1)=0.080000,y(9)=0.950000,f(x,y)=-5.174228132203e-001 x(1)=0.080000,y(10)=1.000000,f(x,y)=-5.036686722854e-001 x(1)=0.080000,y(11)=1.050000,f(x,y)=-4.

38、886112969050e-001 x(1)=0.080000,y(12)=1.100000,f(x,y)=-4.722506870792e-001 x(1)=0.080000,y(13)=1.150000,f(x,y)=-4.545868428078e-001 x(1)=0.080000,y(14)=1.200000,f(x,y)=-4.356197640909e-001 x(1)=0.080000,y(15)=1.250000,f(x,y)=-4.153494509285e-001 x(1)=0.080000,y(16)=1.300000,f(x,y)=-3.937759033206e-0

39、01 x(1)=0.080000,y(17)=1.350000,f(x,y)=-3.708991212673e-001 x(1)=0.080000,y(18)=1.400000,f(x,y)=-3.467191047684e-001 x(1)=0.080000,y(19)=1.450000,f(x,y)=-3.212358538240e-001 x(1)=0.080000,y(20)=1.500000,f(x,y)=-2.944493684341e-001 x(2)=0.160000,y(0)=0.500000,f(x,y)=-5.793772510005e-001 x(2)=0.160000

40、,y(1)=0.550000,f(x,y)=-5.811162261882e-001 x(2)=0.160000,y(2)=0.600000,f(x,y)=-5.814354309237e-001 x(2)=0.160000,y(3)=0.650000,f(x,y)=-5.803348652070e-001 x(2)=0.160000,y(4)=0.700000,f(x,y)=-5.778145290381e-001 x(2)=0.160000,y(5)=0.750000,f(x,y)=-5.738744224169e-001 x(2)=0.160000,y(6)=0.800000,f(x,y

41、)=-5.685145453435e-001 x(2)=0.160000,y(7)=0.850000,f(x,y)=-5.617348978179e-001 x(2)=0.160000,y(8)=0.900000,f(x,y)=-5.535354798400e-001 x(2)=0.160000,y(9)=0.950000,f(x,y)=-5.439162914099e-001 x(2)=0.160000,y(10)=1.000000,f(x,y)=-5.328773325276e-001 x(2)=0.160000,y(11)=1.050000,f(x,y)=-5.204186031931e

42、-001 x(2)=0.160000,y(12)=1.100000,f(x,y)=-5.065401034063e-001 x(2)=0.160000,y(13)=1.150000,f(x,y)=-4.912418331674e-001 x(2)=0.160000,y(14)=1.200000,f(x,y)=-4.745237924762e-001 x(2)=0.160000,y(15)=1.250000,f(x,y)=-4.563859813327e-001 x(2)=0.160000,y(16)=1.300000,f(x,y)=-4.368283997371e-001 x(2)=0.160

43、000,y(17)=1.350000,f(x,y)=-4.158510476892e-001 x(2)=0.160000,y(18)=1.400000,f(x,y)=-3.934539251891e-001 x(2)=0.160000,y(19)=1.450000,f(x,y)=-3.696370322368e-001 x(2)=0.160000,y(20)=1.500000,f(x,y)=-3.444003688322e-001 x(3)=0.240000,y(0)=0.500000,f(x,y)=-5.741362652301e-001 x(3)=0.240000,y(1)=0.55000

44、0,f(x,y)=-5.795336632615e-001 x(3)=0.240000,y(2)=0.600000,f(x,y)=-5.833980237588e-001 x(3)=0.240000,y(3)=0.650000,f(x,y)=-5.857293467218e-001 x(3)=0.240000,y(4)=0.700000,f(x,y)=-5.865276321505e-001 x(3)=0.240000,y(5)=0.750000,f(x,y)=-5.857928800450e-001 x(3)=0.240000,y(6)=0.800000,f(x,y)=-5.83525090

45、4053e-001 x(3)=0.240000,y(7)=0.850000,f(x,y)=-5.797242632314e-001 x(3)=0.240000,y(8)=0.900000,f(x,y)=-5.743903985231e-001 x(3)=0.240000,y(9)=0.950000,f(x,y)=-5.675234962807e-001 x(3)=0.240000,y(10)=1.000000,f(x,y)=-5.591235565040e-001 x(3)=0.240000,y(11)=1.050000,f(x,y)=-5.491905791931e-001 x(3)=0.2

46、40000,y(12)=1.100000,f(x,y)=-5.377245643479e-001 x(3)=0.240000,y(13)=1.150000,f(x,y)=-5.247255119685e-001 x(3)=0.240000,y(14)=1.200000,f(x,y)=-5.101934220548e-001 x(3)=0.240000,y(15)=1.250000,f(x,y)=-4.941282946069e-001 x(3)=0.240000,y(16)=1.300000,f(x,y)=-4.765301296248e-001 x(3)=0.240000,y(17)=1.3

47、50000,f(x,y)=-4.573989271084e-001 x(3)=0.240000,y(18)=1.400000,f(x,y)=-4.367346870578e-001 x(3)=0.240000,y(19)=1.450000,f(x,y)=-4.145374094729e-001 x(3)=0.240000,y(20)=1.500000,f(x,y)=-3.908070943538e-001 x(4)=0.320000,y(0)=0.500000,f(x,y)=-5.668415742753e-001 x(4)=0.320000,y(1)=0.550000,f(x,y)=-5.7

48、57918118812e-001 x(4)=0.320000,y(2)=0.600000,f(x,y)=-5.830990137955e-001 x(4)=0.320000,y(3)=0.650000,f(x,y)=-5.887631800183e-001 x(4)=0.320000,y(4)=0.700000,f(x,y)=-5.927843105496e-001 x(4)=0.320000,y(5)=0.750000,f(x,y)=-5.951624053893e-001 x(4)=0.320000,y(6)=0.800000,f(x,y)=-5.958974645374e-001 x(4

49、)=0.320000,y(7)=0.850000,f(x,y)=-5.949894879940e-001 x(4)=0.320000,y(8)=0.900000,f(x,y)=-5.924384757591e-001 x(4)=0.320000,y(9)=0.950000,f(x,y)=-5.882444278326e-001 x(4)=0.320000,y(10)=1.000000,f(x,y)=-5.824073442145e-001 x(4)=0.320000,y(11)=1.050000,f(x,y)=-5.749272249049e-001 x(4)=0.320000,y(12)=1

50、.100000,f(x,y)=-5.658040699038e-001 x(4)=0.320000,y(13)=1.150000,f(x,y)=-5.550378792111e-001 x(4)=0.320000,y(14)=1.200000,f(x,y)=-5.426286528269e-001 x(4)=0.320000,y(15)=1.250000,f(x,y)=-5.285763907511e-001 x(4)=0.320000,y(16)=1.300000,f(x,y)=-5.128810929837e-001 x(4)=0.320000,y(17)=1.350000,f(x,y)=-4.955427595248e-001 x(4)=0.320000,y(18)=1.400000,f(x,y)=-4.765613903744e-001 x(4)=0.320000,y(19)=1.450000,f(x,y)=-4.559369855324e-001 x(4)=0.320000,y(20)=1.500000,