数值分析计算实习题二

《数值分析》计算实习题二
算法设计方案1．主要计算步骤：
计算函数f(x,y)在拟合所需的节点处的函数值。

将各拟合节点（x i，y j）分别带入非线性方程组
0.5 cos t + u + v + w – x = 2.67
t + 0.5 sin u + v + w – y = 1.07
0.5t + u + cos v + w – x =3.74
t + 0.5u + v + sin w – y =0.79
解非线性方程组得解向量（t ij，u ij，v ij，w ij）。

对数表z(t,u)进行分片二次代数插值，求得对应（t ij，u ij）处的值，即为f(x i，y j) 的值。

对上述拟合节点分别进行x，y最高次数为k（k=0,1,2,3…）次的多项式拟合。

每次拟合后验证误差大小，直到满足要求。

2．求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。

3．对z(t，u)进行插值选择分片二次插值。

4．拟合基函数φr(x)ψs(y)选择为φr(x)=x r，ψs(y)=y s。

拟合系数矩阵c通过连续两次解线性方程组求得。

一．源程序
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
void Doolittle(double *A,int n,int *M)
//功能说明：对n阶矩阵A进行选主元的Doolittle分解
//参数说明：A：欲进行分解的方阵，同时也是返回参数，分解后的结果
// 存储于A中
// n：方阵A的维数
// M；（返回参数）n维向量，记录选主元过程中行交换的次序
{
int i,j,k,t;
double *s;
double Maxs,temp;
s=(double*) calloc(n,sizeof(double));
for(k=0;k<n;k++)
{
for(i=k;i<n;i++)
{
s[i]=A[i*n+k];
for(t=0;t<k;t++) s[i]-= A[i*n+t] * A[t*n+k];
}
Maxs=abs(s[k]); M[k]=k;
for(i=k+1;i<n;i++)
{
if(Maxs<abs(s[i]))
{
Maxs=abs(s[i]);
M[k]=i;
}
}
if(M[k]!=k)
{
for(t=0;t<n;t++)
{
temp=A[k*n+t];
A[k*n+t]=A[M[k]*n+t];
A[M[k]*n+t]=temp;
}
temp=s[k];
s[k]=s[M[k]];
s[M[k]]=temp;
}
if(Maxs<(1e-14))
{
s[k]=1e-14;
printf("%.16e方阵奇异\n",Maxs);
}
A[k*n+k]=s[k];
for(j=k+1;(j<n)&&(k<n-1);j++)
{
for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j];
A[j*n+k]=s[j]/A[k*n+k];
}
}
}
void Solve_LUEquation(double* A,int n,double* b,double* x) //功能说明：解方程LUx＝b，其中L、U共同存储在A中
//参数说明：A：经Doolittle分解后的方阵
// n：方阵A的维数
// b：方程组的右端向量
// x：（返回参数）方程组的解向量
{
int i,t;
for(i=0;i<n;i++)
{
x[i]=b[i];
for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t];
}
for(i=n-1;i>-1;i--)
{
for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t];
x[i]/=A[i*n+i];
}
}
void Transpose(double *A,int m,int n,double* AT)
//功能说明：求m×n阶矩阵A的转置AT
//参数说明：A：已知m×n阶矩阵
// m：A的行数
// n：A的列数
// AT：（返回参数）A的转置矩阵（n×m）
{
int i,j;
for(i=0;i<m;i++)
for(j=0;j<n;j++) AT[j*m+i]=A[i*n+j];
}
void Solve_LEquation(double* A,int n,double* B,double* x,int m) //功能说明：解线性方程组Ax＝B，该函数可对系数矩阵相同
// 而右端向量不同的多个方程组同时求解。

//参数说明：A：方程组系数矩阵
// n：A的维数
// B：m组右端向量构成的n×m矩阵
// x：（返回参数）方程组的解，n×m矩阵
// m：不同右端向量的组数。

{
int* M,i,j;
M=(int*) calloc(n,sizeof(int));
double *BT,*xT,temp;
BT=(double*) calloc(n*m,sizeof(double));
xT=(double*) calloc(n*m,sizeof(double));
//求B的转置BT,使得对应一个方程组的右端系数可以连续存储
//便于函数调用
Transpose(B,n,m,BT);
Doolittle(A,n,M);//将A三角分解
for(i=0;i<m;i++)
{
for(j=0;j<n-1;j++)
{
temp=BT[i*n+j];
BT[i*n+j]=BT[i*n+M[j]];
BT[i*n+M[j]]=temp;
}
//对不同右端向量分别解方程组
Solve_LUEquation(A,n,BT+i*n,xT+i*n);
}
Transpose(xT,m,n,x);
}
void Matrix_Mult(double* A,double* B,int m,int s,int n,double* C) //功能说明：矩阵乘法AB＝C，其中A为m×s阶，B为s×n阶。

//参数说明：A：m×s矩阵
// B：s×n矩阵
// m；A的行数
// s：A的列数，B的行数
// n：B的列数
// C：返回值，m×n矩阵，返回AB的结果
{
int i,j,k;
for(i=0;i<m;i++)
for(j=0;j<n;j++)
{
C[i*n+j]=0;
for(k=0;k<s;k++) C[i*n+j]+=A[i*s+k]*B[k*n+j];
}
}
void Set_NLtoL_JacobiA(double* A,double* x)
//功能说明：求题中非线性方程组对应自变量向量x的雅克比矩阵
//参数说明：A：4×4矩阵
// x：4维向量，非线性方程的自变量（t,u,v,w)
{
A[0]=-0.5*sin(x[0]); A[1]=1;
A[2]=1; A[3]=1;
A[4]=1; A[5]=0.5*cos(x[1]);
A[6]=1; A[7]=1;
A[8]=0.5; A[9]=1;
A[10]=-sin(x[2]); A[11]=1;
A[12]=1; A[13]=0.5;
A[14]=1; A[15]=cos(x[3]);
}
void Set_NLtoL_B(double* B,double* x,double a,double b)
//功能说明：求非线性方程组Newton迭代法的右端式:-F(x)
//参数说明：B：4维向量，对应Newton迭代法的右端式-F(x)
// x：4维向量，非线性方程的自变量（t,u,v,w)
// a：非线性方程的参数，这里对应题中的x
// b：非线性方程的参数，这里对应题中的y
{
B[0]=-(0.5*cos(x[0])+x[1]+x[2]+x[3]-a-2.67);
B[1]=-(x[0]+0.5*sin(x[1])+x[2]+x[3]-b-1.07);
B[2]=-(0.5*x[0]+x[1]+cos(x[2])+x[3]-a-3.74);
B[3]=-(x[0]+0.5*x[1]+x[2]+sin(x[3])-b-0.79);
}
double Vector_FanShu(double *V,int n)
//功能说明：求n维向量V的无穷范数
//返回值：所求范数
{
int i;
double max=0;
for(i=0;i<n;i++)
if(max<abs(V[i])) max=abs(V[i]);
return max;
}
void Solve_NLEquation(double a,double b,double* x)
//功能说明：求解非线性方程组
//参数说明：x：（返回值）4维向量，非线性方程的解向量（t,u,v,w) // a：非线性方程的参数，这里对应题中的x
// b：非线性方程的参数，这里对应题中的y
{
double A[4][4],F[4],detx[4];
double det=1e-12; //求解精度要求
int i,k=0;
int M=5000; //最大迭代次数
//设定迭代初始值
x[0]=0.5; x[1]=1; x[2]=1; x[3]=1;
do
{
Set_NLtoL_B(F,x,a,b);
Set_NLtoL_JacobiA(*A,x);
Solve_LEquation(*A,4,F,detx,1);
if(Vector_FanShu(detx,4)/Vector_FanShu(x,4)<det) return;
for(i=0;i<4;i++) x[i]+=detx[i];
k++;
}while(k<M);
printf("Newton法在该初值不收敛\n");
}
int NearPointIndex(double* v,double a,int n)
//功能说明：查找n维向量V中与常数a最接近的元素的下标
//参数说明：V：n维向量
// n：向量V的维数
//返回值：与a最接近的元素的下标值
{
int i,Index;
double min;
min=abs(v[0]-a);
Index=0;
for(i=1;i<n;i++)
{
if(min>abs(v[i]-a))
{
min=abs(v[i]-a);
Index=i;
}
}
return Index;
}
double ChaZhiZ(double a,double b)
//功能说明：求数表z(t,u)在点（a，b）处的分片二次代数差值
{
double z[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};
double x[6]={ 0, 0.2, 0.4, 0.6, 0.8, 1 };
double y[6]={ 0, 0.4, 0.8, 1.2, 1.6, 2 };
double L[3],L_[3],Z;
int i,j,k,r,t;
i=NearPointIndex(x,a,6);
j=NearPointIndex(y,b,6);
//插值区域边界处插值节点的选取
if(i==0) i=1;
if(i==5) i=4;
if(j==0) j=1;
if(j==5) j=4;
for(k=i-1;k<=i+1;k++)
{
L[k-i+1]=1;
for(t=i-1;t<=i+1;t++)
{
if(t!=k) L[k-i+1]*=((a-x[t])/(x[k]-x[t]));
}
}
for(r=j-1;r<=j+1;r++)
{
L_[r-j+1]=1;
for(t=j-1;t<=j+1;t++)
{
if(t!=r) L_[r-j+1]*=((b-y[t])/(y[r]-y[t]));
}
}
Z=0;
for(k=i-1;k<=i+1;k++)
{
for(r=j-1;r<=j+1;r++)
{
Z+=(L[k-i+1]*L_[r-j+1]*z[k][r]);
}
}
return Z;
}
void NiHe(double*U,double*x,double*y,int m,int n,int p,int q,double*C)
//功能说明：对m×n数表U(x,y)进行二元多项式拟合
//参数说明：U：m×n矩阵，拟合数据点处的函数值
// x；m维向量，第一个自变量的数据节点
// y：n维向量，第二个自变量的数据节点
// m；第一个自变量数据节点个数，也是U的行数
// n：第二个自变量数据节点个数，也是U的列数
// p：构成乘积型基函数中线性无关的x函数的个数，
// 即x最高次数为p－1
// q：构成乘积型基函数中线性无关的y函数的个数，
// 即y最高次数为q－1
// C：p×q矩阵，拟合系数
{
int i,j;
double *B,*BT,*BTB,*G,*GT,*GTG,*BTUG,*temp,*tempT,*CT;
B=(double*) calloc(m*p,sizeof(double));
BT=(double*) calloc(p*m,sizeof(double));
BTB=(double*) calloc(p*p,sizeof(double));
G=(double*) calloc(n*q,sizeof(double));
GT=(double*) calloc(q*n,sizeof(double));
GTG=(double*) calloc(q*q,sizeof(double));
BTUG=(double*) calloc(p*q,sizeof(double));
temp=(double*) calloc(p*n,sizeof(double));
for(i=0;i<m;i++)
{
for(j=0;j<p;j++) B[i*p+j]=pow(x[i],j);
}
for(i=0;i<n;i++)
{
for(j=0;j<q;j++) G[i*q+j]=pow(y[i],j);
}
Transpose(B,m,p,BT);
Transpose(G,n,q,GT);
Matrix_Mult(BT,B,p,m,p,BTB);
Matrix_Mult(GT,G,q,n,q,GTG);
Matrix_Mult(BT,U,p,m,n,temp);
Matrix_Mult(temp,G,p,n,q,BTUG);
free(B);
free(BT);
free(G);
free(GT);
free(temp);
temp=(double*) calloc(p*q,sizeof(double));
Solve_LEquation(BTB,p,BTUG,temp,q);
free(BTB);
free(BTUG);
tempT=(double*) calloc(q*p,sizeof(double));
CT=(double*) calloc(q*p,sizeof(double));
Transpose(temp,p,q,tempT);
Solve_LEquation(GTG,q,tempT,CT,p);
Transpose(CT,q,p,C);
free(temp);
free(CT);
free(tempT);
free(GTG);
}
double Pxy(double *C,int p,int q,double x,double y)
//功能说明：求多项式拟合函数P(x,y)在点（x，y）处的函数值。

//参数说明：C：p×q矩阵，拟合系数矩阵
// p：C的行数
// q：C的列数
{
int i,j;
double sum=0;
for(i=0;i<p;i++)
{
for(j=0;j<q;j++) sum+=(C[i*q+j]*pow(x,i)*pow(y,j));
}
return sum;
}
void main()
{
double NLAns[4]; //非线性方程组变量向量（t,u,v,w）
double f[11][21]; //f(x,y)在拟合节点处的值
double x[11],y[21]; //拟合节点
double *C; //拟合系数矩阵
double det=1e-7; //拟合精度要求
double p,d;
int i,j,k,t;
//设置拟合节点
for(i=0;i<11;i++) x[i]=0.08*i;
for(j=0;j<21;j++) y[j]=0.5+0.05*j;
//求拟合节点处的f(x,y)的值
for(i=0;i<11;i++)
{
for(j=0;j<21;j++)
{
Solve_NLEquation(x[i],y[j],NLAns);
f[i][j]=ChaZhiZ(NLAns[0],NLAns[1]);
}
}
//打印数表f(xi,yi)
printf("数表 f(xi,yi):\n");
for(t=0;t<21/3;t++)
{
printf("┏───┯─────────────────────\────────────┓\n");
printf("│f(x,y)│");
for(i=3*t;i<3*t+3;i++)
{
printf(" y=%4.2f ",y[i]);
}
printf("│\n");
printf("┢───╀─────────────────────\────────────┩\n");
for(j=0;j<11;j++)
{
printf("│x=%4.2f│",x[j]);
for(i=3*t;i<3*t+3;i++) printf("%+.12e ",f[j][i]);
printf("│\n");
}
printf("┗───┻─────────────────────\────────────┚\n");
if(i!=21) printf("\n续表%d:\n",t+1);
}
//k=0,1,2,3....时的拟合
k=0;
do{
C=(double*) calloc((k+1)*(k+1),sizeof(double));
NiHe(*f,x,y,11,21,k+1,k+1,C);
//求在当前k值拟合的节点处误差
d=0;
for(i=0;i<11;i++)
{
for(j=0;j<21;j++)
{
p=Pxy(C,k+1,k+1,x[i],y[j]);
d+=((f[i][j]-p)*(f[i][j]-p));
}
}
printf("\n当k=%d时，σ= %.12e",k,d);
if(d>=det)
free(C);
else
{
printf("<%.0e,满足要求\n",det);
break;
}
}while(++k<11);
//打印拟合系数矩阵c
printf("\n当k=%d时，拟合函数p(x,y)的系数矩阵c：\n",k);
for(t=0;t<(k+1)/3;t++)
{
if(t==0) printf("┏ \n");
if(t==(k+1)/3-1)printf(" \
┓\n");
for(j=0;j<(k+1);j++)
{
if(t==0) printf("│");
for(i=3*t;i<3*t+3;i++) printf("%+.12e ",C[j*(k+1)+i]);
if(t==(k+1)/3-1) printf("│");
printf("\n");
}
if(t==0) printf("┗\n");
if(t==(k+1)/3-1) printf(" \
┚\n\n");
if(i!=k+1) printf("\n续表%d:\n",t+1);
}
//打印(x*,y*)处的拟合逼近情况
printf("观察以下各点的逼近情况\n");
printf("┏─────┮───────────┮─────────\
──┮─────┓\n");
printf("│(x*,y*) │ f(x*,y*) │ p(x*,y*) \ │ |f-p| │\n");
printf("┢─────╋───────────╋─────────\
──╋─────┫\n");
for(i=1;i<=8;i++)
for(j=1;j<=5;j++)
{
Solve_NLEquation(0.1*i,0.5+0.2*j,NLAns);
d=ChaZhiZ(NLAns[0],NLAns[1]);
p=Pxy(C,k+1,k+1,0.1*i,0.5+0.2*j);
printf("│(%3.1f,%3.1f) │ %+.12e │ %+.12e │ \
%f │\n",0.1*i,0.5+0.2*j,d,p,abs(d-p));
}
printf("┗─────┸───────────┻─────────\
──┻─────┙");
getchar();
}
二．计算结果：。