数值分析第三次作业.docx

1、数值分析第三次作业数值分析第三次上机作业算法设计1、求解非线性方程组采用牛顿法解非线性方程组,将题目中给出的当作已知量代入题目给定的非线性方程组，求出与相对应的数组tij，uij 2、分片二次代数插值对所求出的数组tij，uij，通过分片二次代数插值运算，得到与数组t1121，u1121对应的数组z1121，得到二元函数z=，二次插值采用教材给出的分片二次代数插值。3、曲面插值利用x11，y21，z1121建立二维函数表，进行曲面插值计算，逐步提高k值，计算其精度，看其是否满足要求，条件满足则循环结束，并得到曲面拟合的系数矩阵Crs，算法采用教材给出的曲面拟合算法，求出所需矩阵给出，然后按公式

2、进行计算。4、比较逼近的效果观察和逼近的效果时，只需要利用新的点列重复计算二次代数插值，得到与新的插值节点对应的，再与对应的比较即可，求解事可以直接使用（3）中的Crs和k值。源程序(在CodeBlock C/C+集成开发环境下编译通过)#include #include #include #include using namespace std;#define EQUANUM 4#define THTA 1E-12#define SIGMA 1E-7#define ITIME 1000void printMatrix(double *A,int r,int s) for(int i=0;ir

3、;i+) for(int j=0;js;j+) coutAi*s+j ; coutendl; void getFValue(double b,double offset,double x) b1=-(0.5*cos(x1)+x2+x3+x4-offset1); b2=-(x1+0.5*sin(x2)+x3+x4-offset2); b3=-(0.5*x1+x2+cos(x3)+x4-offset3); b4=-(x1+0.5*x2+x3+sin(x4)-offset4);void getFDeri(double coefEQUANUM+1,double x) coef11=-0.5*sin(x

4、1); coef12=1; coef13=1; coef14=1; coef21=1; coef22=0.5*cos(x2); coef23=1; coef24=1; coef31=0.5; coef32=1; coef33=-sin(x3); coef34=1; coef41=1; coef42=0.5; coef43=1; coef44=cos(x4);/解线性方程组void solveEqua(double(*mat)EQUANUM+1,double b,double x) for(int k=1;kEQUANUM;k+) int i=k; for(int j=k;jmatik) i=j

5、; for(int j=k;j=EQUANUM;j+) double temp=matkj; matkj=matij; matij=temp; double temp=bk; bk=bi; bi=temp; for(int i=k+1;i=EQUANUM;i+) double mik=matik/matkk; for(int j=k+1;j=1;k-) xk=0; for(int j=k+1;j=EQUANUM;j+) xk-=matkj*xj; xk+=bk; xk=xk/matkk; double getVecLen(double vec) double temp=0; for(int i

6、=1;i=EQUANUM;i+) if(tempfabs(veci) temp=fabs(veci); return temp;void Newton(double x,double offset) char a; double detxEQUANUM+1; double tempxEQUANUM+1; for(int i=1;i=EQUANUM;i+) xi=1; double FDeriEQUANUM+1EQUANUM+1; double FValueEQUANUM+1; for(int i=1;iITIME;i+) getFDeri(FDeri,x); getFValue(FValue,

7、offset,x); solveEqua(FDeri,FValue,detx); if(getVecLen(detx)/getVecLen(x)=THTA) break; else for(int i=1;i=EQUANUM;i+) xi+=detxi; void biInter(double t21,double u21,double z21,int xs,int ys) double tt6=0,0.2,0.4,0.6,0.8,1; double uu6=0,0.4,0.8,1.2,1.6,2; double zz66=-0.5,-0.34,0.14,0.94,2.06,3.5, -0.4

8、2,-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; double h=0.2,tao=0.4; int n=5,m=5; for(int i=0;ixs;i+) for(int j=0;jys;j+) int xp=0,yp=0; if(tijttn-1-h/2) xp=n-1; else for(int q=2;qttq-h/

9、2)&(tij=ttq+h/2) xp=q; if(uijuum-1-tao/2) yp=n-1; else for(int q=2;quuq-tao/2)&(uij=uuq+tao/2) yp=q; zij=0; for(int k=xp-1;k=xp+1;k+) double lk=1; for(int ti=xp-1;ti=xp+1;ti+) if(ti!=k) lk*=(tij-ttti)/(ttk-ttti); for(int r=yp-1;r=yp+1;r+) double lr=1; for(int ti=yp-1;ti=yp+1;ti+) if(ti!=r) lr*=(uij-

10、uuti)/(uur-uuti); zij+=lk*lr*zzkr; void matrixMut(double *A,double *B,int r,int s,int t,double *result) for(int i=0;ir;i+) for(int j=0;jt;j+) resulti*t+j=0; for(int k=0;ks;k+) resulti*t+j+=Ai*s+k*Bk*t+j; void transpose(double *A,int r,int s,double *result) for(int i=0;ir;i+) for(int j=0;js;j+) resul

11、tj*r+i=Ai*s+j;void copyMatrix(double *A,int r,int s,int wd,double *result) for(int i=0;ir;i+) for(int j=0;js;j+) resulti*s+j=Ai*wd+j;void inversion(double *A,int r,double *result) double copyAr*r; copyMatrix(A,r,r,r,copyA); for(int i=0;ir;i+) for(int j=0;jr;j+) if(i!=j) resulti*r+j=0; else resulti*r

12、+j=1; for(int i=0;ir;i+) int index=i; for(int j=i+1;jfabs(copyAindex*r+i) index=j; if(index!=i) for(int j=0;jr;j+) double temp=copyAindex*r+j; copyAindex*r+j=copyAi*r+j; copyAi*r+j=temp; temp=resultindex*r+j; resultindex*r+j=resulti*r+j; resulti*r+j=temp; double temp=copyAi*r+i; if(temp!=1) for(int

13、j=0;jr;j+) copyAi*r+j/=temp; resulti*r+j/=temp; for(int j=0;jr;j+) if(copyAj*r+i!=0)&(i!=j) temp=copyAj*r+i; for(int k=0;kr;k+) copyAj*r+k-=temp*copyAi*r+k; resultj*r+k-=temp*resulti*r+k; double *surFit(double *z,int &kvalue) int xs=11,ys=21,num=9; double Bxs*num; double Gys*num; double Pxs*ys; doub

14、le TBxs*ys; double TGxs*ys; double BTxs*ys; double GTxs*ys; double BBxs*ys; double *GG=new doublexs*ys; for(int i=0;inum;i+) for(int j=0;jxs;j+) Bj*num+i=pow(0.08*j,i); for(int j=0;jys;j+) Gj*num+i=pow(0.5+0.05*j,i); double sig=0; for(int i=0;inum;i+) sig=0; copyMatrix(B,xs,i+1,num,TB); transpose(TB

15、,xs,i+1,BT); matrixMut(BT,TB,i+1,xs,i+1,BB); inversion(BB,i+1,TB); copyMatrix(G,ys,i+1,num,TG); transpose(TG,ys,i+1,GT); matrixMut(GT,TG,i+1,ys,i+1,GG); inversion(GG,i+1,GT); matrixMut(TB,BT,i+1,i+1,xs,BB); matrixMut(BB,z,i+1,xs,ys,GG); matrixMut(GG,TG,i+1,ys,i+1,BB); matrixMut(BB,GT,i+1,i+1,i+1,GG)

16、; for(int j=0;jxs;j+) for(int k=0;kys;k+) double temp=0; for(int p=0;pi+1;p+) for(int q=0;qi+1;q+) temp+=GGp*(i+1)+q*Bj*num+p*Gk*num+q; Pj*ys+k=temp; sig+=(zj*ys+k-temp)*(zj*ys+k-temp); coutk=isetprecision(11)setiosflags(ios:scientific|ios:uppercase) sigma=sigendl; ofstream out; out.open(shubiao1.tx

17、t,ios:app); outk=isetprecision(11)setiosflags(ios:scientific|ios:uppercase) sigma=sigendl; out.close(); if(sigSIGMA) kvalue=i; out.open(shubiao2.txt,ios:app); outk=isetprecision(11)setiosflags(ios:scientific|ios:uppercase) sigma=sigendl; for(int k=0;k=i;k+) for(int j=0;j=i;j+) outCij=GGk*(i+1)+jendl

18、; coutCij=GGk*(i+1)+jendl; out.close(); return GG; return NULL;int main() double t1121,u1121,z1121,offsetEQUANUM+1,xEQUANUM+1; int kvalue=0; for(int i=0;i=10;i+) for(int j=0;j=20;j+) offset1=2.67+0.08*i; offset2=1.07+0.5+0.05*j; offset3=3.74+0.08*i; offset4=0.79+0.5+0.05*j; Newton(x,offset); tij=x1;

19、 uij=x2; biInter(t,u,z,11,21); ofstream out; out.open(shubiao.txt); for(int i=0;i=10;i+) for(int j=0;j=20;j+) coutresetiosflags(ios:scientific|ios:uppercase)x=0.08*i y=0.5+0.05*j; coutsetprecision(11)setiosflags(ios:scientific|ios:uppercase) f(x,y)=zijendl; outresetiosflags(ios:scientific|ios:upperc

20、ase)x=0.08*i y=0.5+0.05*j; outsetprecision(11)setiosflags(ios:scientific|ios:uppercase) f(x,y)=zijendl; out.close(); double *GG=surFit(z0,kvalue); for(int i=1;i=8;i+) for(int j=1;j=5;j+) offset1=2.67+0.1*i; offset2=1.07+0.5+0.2*j; offset3=3.74+0.1*i; offset4=0.79+0.5+0.2*j; Newton(x,offset); ti-1j-1

21、=x1; ui-1j-1=x2; biInter(t,u,z,8,5); double p85; for(int i=1;i=8;i+) for(int j=1;j=5;j+) double temp=0; for(int ii=0;ii=kvalue;ii+) for(int jj=0;jj=kvalue;jj+) temp+=GGii*(kvalue+1)+jj*pow(0.1*i,ii)*pow(0.5+0.2*j,jj); pi-1j-1=temp; out.open(shubiao3.txt); for(int i=0;i8;i+) for(int j=0;j5;j+) coutre

22、setiosflags(ios:scientific|ios:uppercase); coutxi=0.1*i yj=0.5+0.2*jendl; coutsetprecision(11)setiosflags(ios:scientific|ios:uppercase); coutp(x,y)=pij f(x,y)=zijendl; coutdeta=pij-zijendl; outresetiosflags(ios:scientific|ios:uppercase); outxi=0.1*i yj=0.5+0.2*jendl; outsetprecision(11)setiosflags(i

23、os:scientific|ios:uppercase); outp(x,y)=pij f(x,y)=zijendl; outdeta=pij-zijendl; out.close(); return 0;程序输出1. 数表(xi,yj) 和 f(xi,yj)x=0 y=0.5 f(x,y)=4.46504018481E-001x=0 y=0.55 f(x,y)=3.24683262928E-001x=0 y=0.6 f(x,y)=2.10159686683E-001x=0 y=0.65 f(x,y)=1.03043608316E-001x=0 y=0.7 f(x,y)=3.401895562

24、68E-003x=0 y=0.75 f(x,y)=-8.87358136380E-002x=0 y=0.8 f(x,y)=-1.73371632750E-001x=0 y=0.85 f(x,y)=-2.50534611467E-001x=0 y=0.9 f(x,y)=-3.20276506388E-001x=0 y=0.95 f(x,y)=-3.82668066110E-001x=0 y=1 f(x,y)=-4.37795766738E-001x=0 y=1.05 f(x,y)=-4.85758941444E-001x=0 y=1.1 f(x,y)=-5.26667254884E-001x=0

25、 y=1.15 f(x,y)=-5.60638479797E-001x=0 y=1.2 f(x,y)=-5.87796538768E-001x=0 y=1.25 f(x,y)=-6.08269779090E-001x=0 y=1.3 f(x,y)=-6.22189452876E-001x=0 y=1.35 f(x,y)=-6.29688378186E-001x=0 y=1.4 f(x,y)=-6.30899760003E-001x=0 y=1.45 f(x,y)=-6.25956152545E-001x=0 y=1.5 f(x,y)=-6.14988546609E-001x=0.08 y=0.5 f(x,y)=6.38015226511E-001x=0.08 y=0.55 f(x,y)