《数值分析》第三次大作业
一、算法的设计方案： （一）、总体方案设计：
（1）解非线性方程组。将给定的(xi,yi)当作已知量代入题目给定的非线性方程组，求得与(xi,yi)相对应的数组t[i][j],u[i][j]。
（2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=f(xi,yi)。
（3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k值，并得到曲面拟合的系数矩阵C[r][s]。
（4）观察和p(xi,yi)的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(xi,yi)对应的f(xi,yi)，再与对应的p(xi,yi)比较即可，这里求解
p(xi,yi)可以直接使用（3）中的C[r][s]和k。
（二）具体算法设计：
（1）解非线性方程组
牛顿法解方程组F(x)?0的解x，可采用如下算法： 1）在x附近选取x*(0)*?D，给定精度水平??0和最大迭代次数M。
2）对于k?0,1,?M执行
① 计算F(x(k))和F?(x(k))。 ② 求解关于?x(k)的线性方程组
F?(x(k))?x(k)??F(xk( ))③ 若?x(k)?x(k)?*(k)??，则取x?x，并停止计算；否则转④。
④ 计算x(k?1)?x(k)??x(k)。
⑤ 若k?M，则继续，否则，输出M次迭代不成功的信息，并停止计算。 （2）分片双二次插值
给定已知数表以及需要插值的节点，进行分片二次插值的算法：
??xi?x0?ih(i?0,1,?,n)设已知数表中的点为： ? ，需要插值的节点为(x,y)。
y?y?j?(j?0,1,?,m)?0?j1) 根据(x,y)选择插值节点(xi,yj)：
若x?x1?若y?y1?h2或x?xn?1?h2，插值节点对应取i?1或i?n?1， ，插值节点对应取j?1或i?m?1。
?2或y?yn?1??2hh?x??x?x?,2?i?n?2ii??22 若?
???y??y?y?,2?j?m?2jj??22则选择(xk,yr)(k?i?1,i,i?1;r?j?1,j,j?1)为插值节点。 2）计算
i?1lk(x)??t?i?1t?kj?1x?xtxk?xty?ytyr?ytk(?i?1i,i?,1)
l?r(y)?
r(?j?1j,j?,1)?t?j?1t?r插值多项式的公式为：
i?1j?1 p(x,y)???k?i?1r?j?1lk(x)l?r(y)f(xk,yr)
注：本步进行插值运算的是(t,u)，利用(xi,yj)与(t,u)的对应关系就可以得到z与
(xi,yj)的对应关系。
（3）曲面拟合
根据插值得到的数表xi,yj,f(xi,yj)进行曲面拟合的过程： 1） 根据拟合节点和基底函数写出矩阵B和G：
?(x0)0?0(x)1 B??????(x)0?n(x0)(x1)?(xn)11????1k?(y0)0(x0)??k?0(x1)?(y)1 G????????k??(y)0(xn)??m(y0)(y1)?1???1(ym)1?k(y0)?k?(y1)? ???k?(ym)?2） 计算 C?(BTB)?1BTUG(GTG)?1。
在这里，为了简化计算和编程、避免矩阵求逆，记：
A?(BB)BU，DT?G(GTG)?1
对上面两式进行变形，得到如下两个线性方程组： (BB)A?BU，(GG)D?G
TTTTT?1T通过解上述两个线性方程组，则有：C?AD
kkT3） 对于每一个(xi,yj)， p(xi,yj)?4） 拟合需要达到的精度条件为：
nm**??r?0s?0Crs(xi)(yj)。
rs ????[pi?0j?0(xi,yj)?uij]?102?。
7 其中uij对应着插值得到的数表xi,yj,f(xi,yj)中f(xi,yj)的值。
5） 让k逐步增加，每一次重复执行以上几步，直到
nm*??
??[pi?0j?0(xi,yj)?uij]?102?7 成立。此时的k值就是要求解最小的k。
二、源程序：
#include
#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）*/
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] 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--) { 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; }