?? 數值方法.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 + -