北航数值分析大作业(三)Word文件下载.docx
《北航数值分析大作业(三)Word文件下载.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业(三)Word文件下载.docx（25页珍藏版）》请在冰点文库上搜索。
要求p（x,y）一最小的k值达到以下的精度
其中xi=0.08i，yj=0.5+0.05j。
2．计算f（xi*,yj*）,p（xi*,yj*）（i=1,2,…,8;
j=1,2,…,5）的值，以观察p（x,y）逼近f（x,y）的效果，其中xi*=0.1i,yj*=0.5+0.2j。
二、算法方案
1．使用C++语言实现，使用牛顿迭代法求解非线性方程组，对，，（）的共计11×
21组分别求出非线性方程组的解，即求出与对应的。
均为11×
21的矩阵。
2．由求出的，使用分片二次代数插值法对题中给出的数表进行插值得到。
即得到的11×
21个数值解。
3．k=0时的多项式拟合必然不符合要求，从k=1开始迭代，使用最小二乘法的曲面拟合法对进行拟合，计算在不符合要求的情况下增大。
当时结束计算，输出结果。
4．由3中得到的系数计算的值，再次使用牛顿迭代法对进行求解得到，再次进行二次插值得到结果，以观察逼近的效果。
其中，，。
三、源程序：
#include<
iostream>
vector>
cmath>
algorithm>
iomanip>
#defineN4//方程组未知个数
#defineM6//z,t,u数表阶数
#defineX_Num11
#defineY_Num21//定义数表大小
#defineEPSL1e-12//定义阶数，精度
#defineEPSL21e-7
usingnamespacestd;
typedefvector<
vector<
double>
>
Mat;
//将二维数组简写为Mat
Equation（Matinput）;
//定义求解非线性方程的函数，同时供Inverse，Zxy函数调用
MatInverse（intrank,Matinput2）;
//定义求解逆矩阵的函数
doubleAccuracy（vector<
X_1,vector<
X_2）;
//定义求解近似向量精度的函数
doubleInterpolation（doubleu_1,doublet_1）;
//定义分片代数二次插值函数
MatCrs（vector<
X,vector<
Y,MatU）;
//最小二乘法求解近似表达式系数
MatZxy（vector<
X1,vector<
Y1）;
//定义非线性方程组，调用Equation，Accuracy和Interpolation完成求解
//所有的output应该调整，是否调整为输出到文件为好
voidoutput（vector<
Final1,vector<
Final2,MatFinal3）;
//定义输出函数,输出矩阵
voidoutput2（MatXi）;
doublevector_u[M]={0,0.4,0.8,1.2,1.6,2};
doublevector_t[M]={0,0.2,0.4,0.6,0.8,1};
doublemat_z[M][M]={
{-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},
};
//定义初始数表z,t,u，此处使用数组，而其它矩阵和向量均使用的为vector以及二维vector
voidmain（）
{
inti,j;
vector<
x,y;
x.resize（X_Num）;
y.resize（Y_Num）;
for（i=0;
i<
X_Num;
i++）
{
x[i]=0.08*i;
}
Y_Num;
y[i]=0.5+0.05*i;
}//定义插值节点
MatZ=Zxy（x,y）;
//求出插值点函数值
output（x,y,Z）;
MatCr=Crs（x,y,Z）;
//求出近似多项式
output2（Cr）;
x.resize（8）;
y.resize（5）;
8;
x[i]=0.1*（i+1）;
for（j=0;
j<
5;
j++）
y[j]=0.5+0.2*（j+1）;
}//新插值节点
MatZ2=Zxy（x,y）;
MatP;
P.resize（8）;
P[i].resize（5）;
intk=Cr.size（）;
//利用上一部的Cr获得P值即可
doubletmp;
intm,n;
for（m=0;
m<
m++）
{
for（n=0;
n<
n++）
{
tmp=0;
for（i=0;
k;
{
for（j=0;
{
tmp=tmp+Cr[i][j]*pow（x[m],i）*pow（y[n],j）;
}
}
P[m][n]=tmp;
}
}//使用近似多项式得到的近似值
for（j=0;
{
cout<
<
setprecision
（2）<
scientific<
x[i]<
"
"
;
y[j]<
插值"
setprecision（12）<
Z2[i][j]<
拟合"
P[i][j]<
endl;
}
system（"
pause"
）;
}
Y1）
inti,j,k;
x,y;
x=X1;
y=Y1;
intX_No=x.size（）;
intY_No=y.size（）;
X_1,X_2;
doublewrong;
X_1.resize（N）;
//过渡用于判断误差
X_2.resize（N）;
//此处调试发现x,y值略有差异
MatA;
//使用牛顿法迭代的带求非线性方程
Matt,u,v,w;
MatZ;
//对应t,u的数表Z
A.resize（N）;
N;
A[i].resize（N+1）;
t.resize（X_No）;
u.resize（X_No）;
v.resize（X_No）;
w.resize（X_No）;
Z.resize（X_No）;
X_No;
t[i].resize（Y_No）;
u[i].resize（Y_No）;
v[i].resize（Y_No）;
w[i].resize（Y_No）;
Z[i].resize（Y_No）;
Y_No;
t[i][j]=u[i][j]=w[i][j]=v[i][j]=1;
}//将待求量赋予初值
A[i][j]=1;
A[2][0]=0.5;
A[3][1]=0.5;
i++）//求解对应x,y的t,u,v,w
j=0;
while（j<
Y_No）
A[0][4]=2.67+x[i]-0.5*cos（t[i][j]）-u[i][j]-v[i][j]-w[i][j];
//此处应做修改
A[1][4]=1.07+y[j]-0.5*sin（u[i][j]）-t[i][j]-v[i][j]-w[i][j];
A[2][4]=3.74+x[i]-cos（v[i][j]）-0.5*t[i][j]-u[i][j]-w[i][j];
A[3][4]=0.79+y[j]-sin（w[i][j]）-0.5*u[i][j]-t[i][j]-v[i][j];
A[0][0]=-0.5*sin（t[i][j]）;
A[1][1]=0.5*cos（u[i][j]）;
A[2][2]=-sin（v[i][j]）;
A[3][3]=cos（w[i][j]）;
vector<
Change=Equation（A）;
//调用求解方程得到第一组增量,此处需要注意Change赋值的问题，得到每组增量后怎么处理
for（k=0;
k<
k++）
X_1[k]=Change[k];
X_2[0]=t[i][j];
X_2[1]=u[i][j];
X_2[2]=v[i][j];
X_2[3]=w[i][j];
wrong=Accuracy（X_1,X_2）;
if（wrong<
=EPSL）
t[i][j]=X_2[0];
u[i][j]=X_2[1];
v[i][j]=X_2[2];
w[i][j]=X_2[3];
j++;
else
t[i][j]=X_1[0]+X_2[0];
u[i][j]=X_1[1]+X_2[1];
v[i][j]=X_1[2]+X_2[2];
w[i][j]=X_1[3]+X_2[3];
}
Z[i][j]=Interpolation（u[i][j],t[i][j]）;
//在此处构建循环得到z的*21矩阵，即f（x,y）
returnZ;
Equation（Matinput）//选主元的doolittle分解法，求非线性方程
{
Mata=input;
//获得a
intNum1=a.size（）;
//a的阶数
intNum2=a[0].size（）;
intMax;
vector<
Max2（Num1）;
X（Num1）;
S（Num1）;
doubleQ;
//中转变量
P（Num2）;
inti,j,k,t;
doubletmp;
doubletmp1;
//用于改善矩阵的条件数
for（k=0;
k<
Num1;
k++）
{
Max=k;
for（i=k;
i<
i++）
{
tmp=0;
for（t=0;
t<
k;
t++）
{
tmp=tmp+a[i][t]*a[t][k];
}
S[i]=a[i][k]-tmp;
//计算中间变量s
if（abs（S[i]）>
abs（S[Max]））//应为判断绝对值大小
Max=i;
}//寻找占优的行
P=a[Max];
a[Max]=a[k];
a[k]=P;
//交换行
Q=S[Max];
S[Max]=S[k];
S[k]=Q;
//交换S
a[k][k]=S[k];
for（j=k+1;
j<
j++）
tmp=0;
for（t=0;
k;
t++）
tmp=tmp+a[k][t]*a[t][j];
a[k][j]=a[k][j]-tmp;
}
for（i=k+1;
a[i][k]=S[i]/a[k][k];
}
for（i=1;
doubletmp=0;
for（t=0;
t<
i;
tmp=tmp+a[i][t]*a[t][Num2-1];
a[i][Num2-1]=a[i][Num2-1]-tmp;
}
for（i=0;
Num1;
Max2[i]=abs（a[i][i]）;
for（j=i;
Max2[i]=max（Max2[i],abs（a[i][j]））;
}
Num2;
a[i][j]=a[i][j]/Max2[i];
}//平衡法，用于改善矩阵的条件数
for（i=Num1-1;
i>
=0;
i--）
tmp=0;
for（j=i+1;
tmp=tmp+a[i][j]*X[j];
X[i]=（a[i][Num2-1]-tmp）/a[i][i];
returnX;
X_2）//迭代精度
A=X_1;
B=X_2;
intNum1=A.size（）;
doubletmp1=0;
doubletmp2=0;
tmp1=abs（A[0]）;
for（i=1;
tmp1=max（tmp1,abs（A[i]））;
}//增量的无穷范数
tmp2=abs（B[0]）;
tmp2=max（tmp2,abs（B[i]））;
}//迭代向量的无穷范数
doubleFinal;
Final=tmp1/tmp2;
returnFinal;
doubleInterpolation（doubleu_1,doublet_1）//分片插值函数
doubleu=u_1;
doublet=t_1;
doublez=0;
tmp1;
tmp2;
tmp1.resize（3）;
tmp2.resize（3）;
intk,r;
if（t<
=0.3）
k=0;
if（（0.3<
t）&
&
（t<
=0.5））
k=1;
if（（0.5<
=0.7））
k=2;
if（t>
0.7）
k=3;
}//判断t的范围
if（u<
=0.6）
r=0;
if（（0.6<
u）&
（u<
=1.0））
r=1;
if（（1.0<
=1.4））
r=2;
if（u>
1.4）
r=3;
tmp1[0]=（（t-vector_t[k+1]）*（t-vector_t[k+2]））/（（vector_t[k]-vector_t[k+1]）*（vector_t[k]-vector_t[k+2]））;
tmp1[1]=（（t-vector_t[k]）*（t-vector_t[k+2]））/（（vector_t[k+1]-vector_t[k]）*（vector_t[k+1]-vector_t[k+2]））;
tmp1[2]=（（t-vector_t[k]）*（t-vector_t[k+1]））/（（vector_t[k+2]-vector_t[k+1]）*（vector_t[k+2]-vector_t[k]））;
tmp2[0]=（（u-vector_u[r+1]）*（u-vector_u[r+2]））/（（vector_u[r]-vector_u[r+1]）*（vector_u[r]-vector_u[r+2]））;
tmp2[1]=（（u-vector_u[r]）*（u-vector_u[r+2]））/（（vector_u[r+1]-vector_u[r]）*（vector_u[r+1]-vector_u[r+2]））;
tmp2[2]=（（u-vector_u[r]）*（u-vector_u[r+1]））/（（vector_u[r+2]-vector_u[r+1]）*（vector_u[r+2]-vector_u[r]））;
3;
tmp=1;
tmp=tmp*mat_z[k+i][r+j]*tmp1[i]*tmp2[j];
z=z+tmp;
returnz;
MatInverse（intrank,Matinput2）//用于矩阵求逆，并返回逆矩阵
MatCopy=input2;
intR=rank;
MatX;
X.resize（R）;
R;
X[i].resize（R）;
for（k=0;
if（k==i）
Copy[k].push_back
（1）;
Copy[k].push_back（0）;
}//该法求解逆矩阵的原理，为在input2方阵后，添加列向量，全体列向量组成单位阵，则各解向量组成逆矩阵
vector<
x=Equation（Copy）;
X[j][i]=x[j];
Copy[k].pop_back（）;
returnX;
Y,MatU）//此处可能重定义,XY为x,y的坐标
x=X;
y=Y;
intX_No=x.size（）;
intY_No=y.size（）;
inti,j,k,l,m,n;
MatB,G;
B.resize（X_No）;
G.resize（Y_No）;
MatU_0=U;
for（k=1;
k++）//增加迭代次数的限制，并设置误差的分析
MatC,C_1;
MatU_1;
//U的过渡矩阵,B的转制称U_0,阶数为k+1乘Y_No
U_1.resize（k+1）;
MatB_1,G_1;
C.resize（k+1）;
C_1.resize（k+1）;
B_1.resize（k+1）;
G_1.resize（k+1）;
for（i=0;
k+1;
C[i].resize（k+1）;
C_1[i].resize（k+1）;
U_1[i].
升级会员 