📄 数值方法.cpp
字号:
// 数值方法.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include "iostream.h"
#include "math.h"
inline double& Data(double*A,int n, int x, int y)
{//获得矩阵的x行j列，完全的矩阵
return *(A+x*n+y);
}
inline double& DataSym(double *A,int n, int x, int y)
{//获得矩阵的x行j列，对称的矩阵
if(x>y)
{
int t;
t=x;x=y;y=t;
}
int temp=0;
for(int i=1;i<x;i++)
temp+=n--;
temp+=(y-i);
return *(A+temp);
}
//n拉格朗日插值计算公式Ln(x)
double Ln( double xx , int n , double* x , double* y );
//n次牛顿向前插值计算公式Nnbefore(x)
double Nnbefore( double xx, int n , double x0 , double h , double* y );
//n次牛顿向前插值计算公式Nnbehind(x)
double Nnbehind( double xx, int n , double x0 , double h , double* y );
//杜力特尔三角分解求解线性方程组
void DLTR(int n , double* A , const double* b, double* x);
//乔立斯基三角分解求解对称正定方程组
void QLSJ(int n , double* A , const double* b, double* x);
//输入数据函数，将数据输入到指针A所指向的空间
void InputDataToMatrixA( int n , double *A);
//输入向量函数，将数据输入到指针b所指向的空间
void InputDataVector(int n , double* b);
//输出矩阵A中的数据按照n维n列形式输出
void OutputDataMatrix(int n, const double* A);
//输出向量函数
void OutputDataVector(int n ,const double* x);
//输入对称矩阵函数，将数据输入到指针A所指向的空间
void InputDataToSymMatrixA( int n , double *A);
//输出对称矩阵A中的数据按照n维n列形式输出
void OutputDataSymMatrix(int n, double* A);
//高斯赛德尔迭代法
void GSSDE(double *a,double* b,int n,double * x,double eps);
//雅可比迭代法
void YKB(double *a,double* b,int n,double * x,double eps);
int main(int argc, char* argv[])
{
cout.setf(ios::fixed,ios::floatfield);
cout.width(10);
//*---------------------------------------------------*/
cout<<"/*---------------------------------------------------*/";
//测试拉格朗日插值的正确性
double aa[3]={0.2,0.3,0.4},bb[3]={1.2214,1.3499,1.4918};
cout<<"67页第一大题(1)"<<endl;
cout<<"为了验证正确性，用67页2题作为检验数据"<<endl;
cout<<"其中f(0.2)=1.2214;f(0.3)=1.3499;f(0.4)=1.4918。"<<endl;
cout<<"通过计算f(0.33)拉格朗日插值的结果是"<<Ln(0.33,2,aa,bb)<<endl;
//*---------------------------------------------------*/
cout<<"/*---------------------------------------------------*/";
//验证牛顿向前插值的正确性
double cc[3]={0.47943,0.56464,0.64422};
cout<<"67页第一大题(2)"<<endl;
cout<<"为了验证正确性，用36页l例题6作为检验数据"<<endl;
cout<<"其中f(0.5)=0.47943;f(0.6)=0.56464;f(0.7)=0.64422。"<<endl;
cout<<"通过计算f(0.57891)牛顿向前插值的结果是"<<Nnbefore(0.57891,2,0.5,0.1,cc)<<endl;
//*---------------------------------------------------*/
cout<<"/*---------------------------------------------------*/";
//验证牛顿向后插值的正确性
double dd[3]={0.38942,0.47943,0.56464};
cout<<"67页第一大题(3)"<<endl;
cout<<"为了验证正确性，用36页l例题6作为检验数据"<<endl;
cout<<"其中f(0.4)=0.38942;f(0.5)=0.4793;f(0.6)=0.56464，\n步长是0.1，初始xo=0.4，n=2次插值"<<endl;
cout<<"通过计算f(0.57891)牛顿向后插值的结果是"<<Nnbehind(0.57891,2,0.4,0.1,dd)<<endl;
//*---------------------------------------------------*/
cout<<"/*---------------------------------------------------*/";
//67页的计算的第二题(1)
cout<<"67页第二大题1"<<endl;
cout<<"下面利用编制的(3)字程序计算ln1.54和ln1.98"<<endl;
cout<<"其中步长0.1，计算ln1.54的初始值是1.5，n=2次插值\n计算ln1.98的初始值是1.8"<<endl;
double qq[3];
for(int i=0;i<3;i++)
qq[i]=log(1.5+0.1*i);
cout<<"通过计算ln1.54的牛顿向后插值的结果是"<<Nnbehind(1.54,2,1.5,0.1,qq)<<endl<<"数学表中的结果是"<<log(1.54)<<endl;
for(int g=0;g<3;g++)
qq[g]=log(1.8+0.1*g);
cout<<"通过计算ln1.98的牛顿向后插值的结果是"<<Nnbehind(1.98,2,1.8,0.1,qq)<<endl<<"数学表中的结果是"<<log(1.98)<<endl;
//*---------------------------------------------------*/
//200页第三大题
cout<<"/*---------------------------------------------------*/";
cout<<"199页的第一大题的(3)\n将不用具体数据来检验，将通过三、五题来测试"<<endl;
cout<<"/*---------------------------------------------------*/";
cout<<"200页第三大题"<<endl;
double s[4][4],bbb[4],xxx[4]={0,0,0,0};
InputDataToSymMatrixA(4,s[0]);
cout<<"输入向量b："<<endl;
InputDataVector(4,bbb);
QLSJ(4,s[0],bbb,xxx);
cout<<"利用平方根法解得的解，如下："<<endl;
OutputDataVector(4,xxx);
//*---------------------------------------------------*/
cout<<"/*---------------------------------------------------*/";
cout<<"200页第五大题"<<endl;
cout<<"取b1={1,1,1,1,1,1}T";
double pp[6][6]={0};
InputDataToMatrixA(6,pp[0]);
double b[6]={1,1,1,1,1,1},x[6];
DLTR(6,pp[0],b,x);
cout<<"经过杜力特尔三角分解后矩阵L如下："<<endl;
for( i=0;i<6;i++)
{
for(int j=0;j<6;j++)
if(i==j)
{
cout.width(10);
cout<<(float)1.00;
}
else
if(i>j)
{
cout.width(10);
cout<<Data(pp[0],6,i,j);
}
else
{
cout.width(10);
cout<<(float)0.0;
}
cout<<endl;
}
cout<<"经过杜力特尔三角分解后矩阵U如下："<<endl;
for( i=0;i<6;i++)
{
for(int j=0;j<6;j++)
if(i<=j)
{
cout.width(10);
cout<<Data(pp[0],6,i,j);
}
else
{
cout.width(10);
cout<<(float)0.0;
}
cout<<endl;
}
cout<<"经过计算求解b1，x1如下(输出的是它们的转置)："<<endl;
OutputDataVector(6,b);
OutputDataVector(6,x);
for(int it=2;it<=10;it++)
{
double ss[6][6]={{1,2,4,7,11,16},{2,3,5,8,12,17}
,{4,5,6,9,13,18},{7,8,9,1014,19}
,{11,12,13,14,15,20},{16,17,8,19,20,21}};
double temp=fabs(x[0]);
for(int i=1;i<6;i++)
if(fabs(x[i])>temp)
temp=fabs(x[i]);
for(i=0;i<6;i++)
b[i]=x[i]/temp;
DLTR(6,ss[0],b,x);
cout<<"经过计算求解b"<<it<<"，x"<<it<<"如下(输出的是它们的转置)："<<endl;
OutputDataVector(6,b);
OutputDataVector(6,x);
}
//*---------------------------------------------------*/
cout<<"//*---------------------------------------------------*//";
cout<<"253页第一大题\n将不用具体数据来检验，将通过二题(1)、(2)来测试"<<endl;
cout<<"//*---------------------------------------------------*//";
cout<<"253页第二题(1)、(2)题"<<endl;
InputDataToMatrixA(6,pp[0]);
cout<<"输入向量b"<<endl;
InputDataVector(6,b);
cout<<"输入初始解"<<endl;
InputDataVector(6,x);
YKB(pp[0],b,6,x,0.00001);
cout<<"通过雅可比迭代法计算的解如下："<<endl;
OutputDataVector(6,x);
x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=1;
GSSDE(pp[0],b,6,x,0.00001);
cout<<"通过高斯赛德尔迭代法计算的解如下："<<endl;
OutputDataVector(6,x);
cout<<"//*---------------------------------------------------*//";
getchar();
return 0;
}
double Ln( double xx , int n , double* x , double* y )
{ //n拉格朗日插值计算公式Ln(x)
//并不进行安全检查不管x,y的指针所指的地方是否存在数据
double result=0;
for(int k=0 ; k<=n ;k++)
{
double temp=1;
for( int i=0 ; i<=n ;i++)
if(k!=i)
temp*=(xx-x[i])/(x[k]-x[i]);
result+=temp*y[k];
}
return result;
}
void CreatTableBefore(const double* y , int n , double result[])
{//在result[1]存储一阶向前插分
for(int i=1 ; i<=n ; i ++ )
result[i] = y[i]-y[i-1];
for( int j=n-1 ; j>=1 ;j--)
for( int k=j; k>=1 ; k--)
result[k+1]-=result[k];
}
double Nnbefore( double xx, int n , double x0 , double h , double* y )
{ //n次牛顿向前插值计算公式Nnbefore(x)
//创建向前差分表
double* table =new double[n+1];
CreatTableBefore(y,n,table);
double t=(xx-x0)/h,result=y[0];
for(int i=1;i<=n;i++)
{
double temp=1;
for(int j=0;j<i;j++)
temp*=(t-j)/(j+1);
result+=temp*table[i];
}
return result;
}
double Nnbehind( double xx, int n , double x0 , double h , double* y )
{ //n次牛顿向前插值计算公式Nnbehind(x)
x0+=n*h;
for(int i=0;i<=(n+1)/2-1;i++)
{
double t;
t=y[i];
y[i]=y[n-i];
y[n-i]=t;
}
h*=-1;
return Nnbefore( xx, n , x0 , h , y );
}
void InputDataToMatrixA( int n , double* A)
{ //参数中的n代表了要输入的矩阵的维数
cout<<"请输入n维n列的矩阵的元素(按照从上行到下行从左列到右列):"<<endl;
if( A==0)
return;
for(int i=0;i<n*n;i++)
cin>>A[i];
cout<<'\n'<<"你输入的矩阵如下所示:"<<endl;
for(i=0;i<n*n;i++)
{
if(i%n==0)
cout<<endl;
cout.width(10);
cout<<A[i];
}
cout<<endl;
}
void InputDataVector(int n , double* b)
{
// cout<<"请输入n维列向量b："<<endl;
for(int i=0;i<n;i++)
cin>>b[i];
cout<<endl;
}
void OutputDataVector(int n ,const double* x)
{//n表示向量的维数，b表示要输出的向量结果的地址
// cout<<"计算结果向量如下："<<endl;
for(int i=0;i<n;i++)
{
cout.width(14);
cout<<x[i];
}
cout<<endl;
}
void OutputDataMatrix(int n, const double* A)
{
cout<<"计算结果的矩阵如下："<<endl;
for(int i=0;i<n*n;i++)
{
if(i%n==0)cout<<endl;
cout.width(10);
cout<<A[i];
}
}
void InputDataToSymMatrixA( int n , double *A)
{
cout<<"请输入n维n列的矩阵的元素(按照从上行到下行从左列到右列，只需要输入上三角元素即可):"<<endl;
if( A==0)
return;
for(int i=0;i<(n*n+n)/2;i++)
cin>>A[i];
cout<<'\n'<<"你输入的矩阵如下所示:"<<endl;
for(i=0;i<n*n;i++)
{
if(i%n==0)
cout<<endl;
cout.width(10);
cout<<DataSym(A,n,(i+1)/n+((i+1)%n==0?0:1),((i+1)%n==0?4:(i+1)%n))<<" ";
}
cout<<endl;
}
void OutputDataSymMatrix(int n, double* A)
{
for(int i=0;i<n*n;i++)
{
if(i%n==0)
cout<<endl;
cout.width(10);
cout<<DataSym(A,n,(i+1)/n+((i+1)%n==0?0:1),((i+1)%n==0?4:(i+1)%n))<<" ";
}
cout<<endl;
}
void DLTR(int n , double* A , const double* b, double* x)
{ //杜力特尔三角分解求解线性方程组
//首先利用三角分解的紧凑格式分解矩阵
for( int i=2;i<=n;i++)
Data(A,n,i-1,0)=Data(A,n,i-1,0)/Data(A,n,0,0);
for( int r = 2;r<=n;r++)
{
for( int i=r;i<=n;i++)
{
double temp=0;
for(int k=1;k<=r-1;k++)
temp+=Data(A,n,k-1,i-1)*(r==k?1:Data(A,n,r-1,k-1));
Data(A,n,r-1,i-1)-=temp;
}
for( i=r+1;i<=n;i++)
{
double temp=0;
for(int k=1;k<=r-1;k++)
temp+=Data(A,n,k-1,r-1)*(i==k?1:Data(A,n,i-1,k-1));
Data(A,n,i-1,r-1)-=temp;
Data(A,n,i-1,r-1)/=Data(A,n,r-1,r-1);
}
}
//开始解方程组Ly=b
x[1-1]=b[1-1];
for( i=2;i<=n;i++)
{
double temp=0;
for(int k=1;k<=i-1;k++)
temp+=(i==k?1:Data(A,n,i-1,k-1))*x[k-1];
x[i-1]=b[i-1]-temp;
}
//开始解Ux=y
x[n-1]=x[n-1]/Data(A,n,n-1,n-1);
for(i=n-1;i>=1;i--)
{
double temp=0;
for(int k=i+1;k<=n;k++)
temp+=Data(A,n,i-1,k-1)*x[k-1];
x[i-1]=(x[i-1]-temp)/Data(A,n,i-1,i-1);
}
}
void QLSJ(int n , double* A , const double* b, double* x)
{
//生成矩阵L飘
for(int j=1;j<=n;j++)
{
double temp=0;
for(int k=1;k<=j-1;k++)
temp+=DataSym(A,n,j,k)*DataSym(A,n,j,k);
DataSym(A,n,j,j)=sqrt(DataSym(A,n,j,j)-temp);
for(int i=j+1;i<=n;i++)
{
temp=0;
for(k=1;k<=j-1;k++)
temp+=DataSym(A,n,j,k)*DataSym(A,n,i,k);
DataSym(A,n,i,j)=(DataSym(A,n,i,j)-temp)/DataSym(A,n,j,j);
}
}
//Ly=b
x[1-1]=b[1-1]/DataSym(A,n,1,1);
for(int i=2;i<=n;i++)
{
double temp=0;
for(int k=1;k<=i-1;k++)
temp+=DataSym(A,n,i,k)*x[k-1];
x[i-1]=(b[i-1]-temp)/DataSym(A,n,i,i);
}
x[n-1]=x[n-1]/DataSym(A,n,n,n);
for(i=n-1;i>=1;i--)
{
double temp=0;
for(int k=i+1;k<=n;k++)
temp+=DataSym(A,n,k,i)*x[k-1];
x[i-1]=(x[i-1]-temp)/DataSym(A,n,i,i);
}
}
void GSSDE(double *a,double* b,int n,double * x,double eps)
//高斯赛德尔迭代法
{
double fan2=(eps+1)*(eps+1);
int nn=0;
while(sqrt(fan2)>eps)
{
fan2=0;
for(int i=1;i<=n;i++)
{
double xx=x[i-1];
double temp=0;
for(int j=1;j<=n;j++)
if(j!=i)
temp+=Data(a,n,i-1,j-1)*x[j-1];
x[i-1]=(b[i-1]-temp)/Data(a,n,i-1,i-1);
fan2+=(x[i-1]-xx)*(x[i-1]-xx);
}
nn++;
}
cout<<"一共迭代了"<<n<<"次";
}
void YKB(double *a,double* b,int n,double * x,double eps)
//雅可比迭代法
{
double* p=new double[n];
double fan2=(eps+1)*(eps+1);
int nn=0;
while(sqrt(fan2)>eps)
{
for(int tt=0;tt<n;tt++)
p[tt]=x[tt];
fan2=0;
for(int i=1;i<=n;i++)
{
double temp=0;
for(int j=1;j<=n;j++)
if(j!=i)
temp+=Data(a,n,i-1,j-1)*p[j-1];
x[i-1]=(b[i-1]-temp)/Data(a,n,i-1,i-1);
fan2+=(x[i-1]-p[i-1])*(x[i-1]-p[i-1]);
}
nn++;
}
delete [] p;
cout<<"一共迭代了"<<n<<"次";
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -