{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 t0[6]={0,0.2,0.4,0.6,0.8,1.0}; double u0[6]={0,0.4,0.8,1.2,1.6,2.0}; double temp1,temp2; for(i=0;i<=10;i++) /*选取插值节点*/ { for(j=0;j<=20;j++) { if(t[i][j]<=0.3) s[i][j]=1; else if(t[i][j]>0.7) s[i][j]=4; else { for(m=2;m<=3;m++) { if((t[i][j]>0.2*m-0.1)&&(t[i][j]<=0.2*m+0.1)) { s[i][j]=m; } } } } } for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { f(u[i][j]>1.4) r[i][j]=4; else { for(m=2;m<=3;m++) { if((u[i][j]>0.4*m-0.2)&&(u[i][j]<=0.4*m+0.2)) { r[i][j]=m; } } } } } for(i=0;i<=10;i++) /*插值运算*/ { for(j=0;j<=20;j++) { z[i][j]=0; for(i1=s[i][j]-1;i1<=s[i][j]+1;i1++)
dF[0][0]=-0.5*sin(X[0]); /*求解F'（x）*/ dF[0][1]=1; dF[0][2]=1; dF[0][3]=1; dF[1][0]=1; dF[1][1]=0.5*cos(X[1]); dF[1][2]=1; dF[1][3]=1; dF[2][0]=0.5; dF[2][1]=1; dF[2][2]=-sin(X[2]); dF[2][3]=1; dF[3][0]=1; dF[3][1]=0.5; dF[3][2]=1; dF[3][3]=cos(X[3]); /*高斯选主元消去法求解Δx*/ for(k=0;k<3;k++) { ik=k; for(l=k;l<=3;l++) {if(dF[ik][k]<dF[l][k]) ik=l; } /*选主元*/ temp=0; temp=F[ik]; F[ik]=F[k]; F[k]=temp; for(l=k;l<=3;l++) { temp=0; temp=dF[ik][l]; dF[ik][l]=dF[k][l]; dF[k][l]=temp; } for(l=k+1;l<=3;l++) { m=dF[l][k]/dF[k][k]; F[l]=F[l]-m*F[k]; for(p=k+1;p<=3;p++) {dF[l][p]=dF[l][p]-m*dF[k][p];} } /*消去过程*/ } dx[3]=-F[3]/dF[3][3]; for(k=2;k>=0;k--)
{ for(j1=r[i][j]-1;j1<=r[i][j]+1;j1++) { temp1=1.0; for(m=s[i][j]-1;m<=s[i][j]+1;m++) { if(m!=i1) temp1*=(t[i][j]-t0[m])/(t0[i1]-t0[m]); } temp2=1.0; for(m=r[i][j]-1;m<=r[i][j]+1;m++) { if(m!=j1) temp2*=(u[i][j]-u0[m])/(u0[j1]-u0[m]); } z[i][j]+=temp1*temp2*z0[i1][j1]; } } } } } void qmnh() /*曲面拟合子程序*/ { int i,j,k,m,l,i1,j1,ik; double A[N][21],D[N][21],B[11][N],Bt[N][11],BtB[N][N],BtU[N][21],BtB1[N][N]; double G[21][N],Gt[N][21],GtG[N][N],GtG1[N][N]; double sigma,p[11][21],temp,q; printf("选择过程的k和sigma值:\n"); k=0; sigma=1; /*选取初值，使循环开始*/ /*求解系数矩阵Crs*/ while(sigma>1e-7) { for(i=0;i<=10;i++) { for(j=0;j<=k;j++) { B[i][j]=pow(x[i],j); Bt[j][i]=B[i][j]; } } for(i=0;i<=k;i++) { for(j=0;j<=k;j++) { temp=0; for(l=0;l<=10;l++) { temp+=Bt[i][l]*B[l][j]; } BtB[i][j]=temp; } } for(i=0;i<=k;i++) { for(j=0;j<=20;j++)
《数值分析》第三次大作业
1、 算法的设计方案： （一）、总体方案设计：
（1）解非线性方程组。将给定的当作已知量代入题目给定的非线性方程组，求得与相对应 的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]] 对应的数组z[11][21]，得到二元函数z=。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的 插值节点对应的，再与对应的比较即可，这里求解可以直接使用（3）中的C[r][s]和k。
{ temp=0; for(l=k+1;l<=3;l++) {temp=temp+dF[k][l]*dx[l]/dF[k][k];} dx[k]=-F[k]/dF[k][k]-temp; } temp=0; for(l=0;l<=3;l++) /*求解矩阵范数，用无穷范数*/ { if(temp<fabs(dx[l])) temp=fabs(dx[l]); } fx=temp; temp=0; for(l=0;l<=3;l++) { if(temp<fabs(X[l])) temp=fabs(X[l]); } fX=temp; if(fabs(fx/fX)<Epsilon1) /*判断是否成立*/ { t[i][j]=X[0]; u[i][j]=X[1]; goto loop4;} else { for(l=0;l<=3;l++) {X[l]=X[l]+dx[l];} n=n+1; goto loop3;} } loop3:{if(n<M) /*判断是否超出规定迭代次数*/ goto loop1; else printf("迭代不成功\n"); goto loop4; } loop4:{continue;} } } } void fpeccz(double t[11][21],double u[11][21])/*分片二次代数插值子程序*/ { int s[11][21],r[11][21]; int i,j,i1,j1,m; double z0[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},
对上面两式进行变形，得到如下两个线性方程组： ， 通过解上述两个线性方程组，则有： 3） 对于每一个， 。 4） 拟合需要达到的精度条件为： 。 其中对应着插值得到的数表中的值。 5） 让k逐步增加，每一次重复执行以上几步，直到 成立。此时的k值就是要求解最小的k。
2、 源程序：
#include<stdio.h> #include<iostream> #include<stdlib.h> #include<math.h> #include<float.h> #include<iomanip> #define Epsilon1 1e-12 /*解线性方程组时近似解向量的精度*/ #define M 200 /*解线性方程组时的最大迭代次数*/ #define N 10 /*求解迭代次数时假设的k的最大值，用于定义包含k的存储空间*/ void Newton(); /*牛顿法求解非线性方程组子程序*/ void fpeccz(); /*分片二次代数插值子程序*/ void qmnh(); /*曲面拟合子程序*/ void duibi(); /*对比 和p逼近效果的子程序*/ double x[11],y[21],t[11][21],u[11][21];/*定义全局变量*/ double z[11][21],C[10][10]; double kz; void Newton(double x[11],double y[21])/*牛顿法求解非线性方程组子程序*/ { double X[4],dx[4],F[4],dF[4][4],temp,m,fx,fX; int i,j,k,l,p,ik,n; for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { X[0]=1; /*选取迭代初始向量，四个分别代表t，u，v，w*/ X[1]=1; X[2]=1; X[3]=1; n=0; loop1:{ F[0]=0.5*cos(X[0])+X[1]+X[2]+X[3]-x[i]-2.67; F[1]=X[0]+0.5*sin(X[1])+X[2]+X[3]-y[j]-1.07; F[2]=0.5*X[0]+X[1]+cos(X[2])+X[3]-x[i]-3.74; F[3]=X[0]+0.5*X[1]+X[2]+sin(X[3])-y[j]-0.79; /*求解F（x）*/