时间序列分析

所谓的直接法，就是经过有限此的运算即可求得方程组的精确解的方法（如果没有舍入误差的话）。他的优点是可以直接得到方程组得精确解，适合与系数矩阵阶数不太高得时候，但是程序较复杂。

高斯消元法的步骤分为两个过程：消元和回代。这个方法在解线性方程组中是常用的一种方法。但是，当它用于程序设计中时可能会存在一个不稳定、不可靠的问题，不能对所有的系数矩阵都能适用。它是依次地对系数矩阵的对角变量进行消元。这就要求这个系数矩阵必须是严格对角占优的。只有这样才能保证消元过程中高斯消元的主元不等0。

如果，在程序设计中直接用高斯消元法将极可能发生除零中断。而且，算法无法依次消元。

要实现一个计算稳定，能应对所有“严格对角占优”和“非严格对角占优”的系数矩阵，并且计算灵巧、高效，必须对高斯消元法进行一些改进。

雅克比迭代（Jacobi）、高斯－赛德尔迭代（Gauss－Seidel）具体程序实现如下：

//雅克比迭代

#include <iostream>

#include <string>

#include <vector>

#include <iomanip>

using namespace std;

const int n=3;????????? //设置方程组的维

float a[n][n],x[n],b[n];

void input_data() //输入方程组的相关数据

{??

??? cout<<"输入方程组的系数矩阵a["<<n<<"]["<<n<<"]:"<<endl;

??? for(int i=0;i<n;i++)

??????? for(int j=0;j<n;j++)

??????????? cin>>a[i][j];

??? cout<<"输入x[1,...,"<<n<<"]的初值:"<<endl;

??? for(i=0;i<n;i++)???????

??????????? cin>>x[i];

??? cout<<"输入b[1,...,"<<n<<"]的初值:"<<endl;

??? for(i=0;i<n;i++)???????

??????????? cin>>b[i];?

}

void output_equation()?//输出方程组

{

??? for(int i=0;i<n;i++)

??? {

??????? for(int j=0;j<n;j++)

??????? {

??????????? if(j>0&&a[i][j]>0)

??????????????? cout<<"+";

??????????? cout<<a[i][j]<<"?x["<<j+1<<"]?";

??????? }

??????? cout<<"="<<b[i]<<endl;

??? }

}

void main ()

{

??? input_data();

??? int counter;??????????? //设置计算时迭代的次数

??? cout<<"输入迭代的次数"<<endl;

??? cin>>counter;

??? system("cls");

??? output_equation();

??? int k=0,i=0;??????????? //k用来统计迭代的次数;i用来表示x的序列

??? float m_num1;???? //m_x用来保存x的原先值;m_num1用来计算j!=i的数据之和

??? cout<<endl;

??? cout<<"k\t";

??? for(int j=0;j<n;j++)

??????? cout<<"x["<<j+1<<"]\t\t\t";

??? cout<<endl;

??? cout<<k<<"\t";

??? k++;

??? cout << setiosflags(ios::fixed);?? //和下一句结合起来用来设定小数位的位数

??? cout << setprecision(10);

??? float x_next[n];

??? for(j=0;j<n;j++)

??? {

??????? cout<<x[j]<<"\t\t";

??? }

??? cout<<endl;

??? do

??? {??????

??????? m_num1=0;

??????? for(j=0;j<n;j++)

??????? {

??????????? if(j!=i)

??????????????? m_num1=m_num1+a[i][j]*x[j];????

??????? }

??????? x_next[i]=(b[i]-m_num1)/a[i][i];???????

??????? if(i<n-1)

??????????? i++;

??????? else

??????? {

??????????? for(j=0;j<n;j++)

??????????????? x[j]=x_next[j];

??????????? cout<<k<<"\t";

??????????? for(j=0;j<n;j++)

??????????????? cout<<x[j]<<"\t\t";

??????????? cout<<endl;

??????????? i=0;

??????????? k++;

??????? };?????

??? }while(k<=counter);

}

//高斯－赛德尔迭代

#include <iostream>
#include <cmath>
using namespace std;
int a,b,m;
double *x0;
void Jacobi(double **c,double *d,int n,double eps);
void Gauss(double **c,double *d,int n,double eps);
void main()
{
int n;
double **A,*B;
double e;
cout<<"请选择求方程组的迭代方法!雅可比选0，高斯-赛德尔选1!"<<endl;
cin>>n;
cout<<"输入方程组的变量的个数以及方程的个数!"<<endl;
cin>>a>>b;
A=new double*[b];
for(int i=0;i<b;i++)
{
A[i]=new double[a];
}
B=new double[b];
x0=new double[a];
cout<<"输入每个方程组的变量的系数以及方程右端的值!"<<endl;
for(int k=0;k<b;k++)
{
for(int j=0;j<a;j++)
{
cin>>A[k][j];
cout<<"A["<<k<<j<<"]="<<A[k][j]<<endl;
}
cin>>B[k];
cout<<"B["<<k<<"]="<<B[k]<<endl;
}
cout<<"输入方程组迭代的次数及所要求的精度!"<<endl;
cin>>m>>e;
cout<<"输入方程组迭代的初值!"<<endl;
for(int j=0;j<a;j++)
{
cin>>x0[j];
}
switch (n)
{
case 0:Jacobi(A,B,m,e);
break;
case 1:Gauss(A,B,m,e);
break;
default:cout<<"你没有选择求解方程组的一种方法!!"<<endl;
break;
}
}

void Jacobi(double **c,double *d,int n,double eps)
{
int k,i;
double *y = new double[a],*x=new double[a],s,temp=0.0;
k=1;
while(1)
{
temp = 0.0;
for(i=0;i<a;i++)
{
s=0.0;
for(int j=0;j<a;j++)
{
if(j!=i)
{
s+=c[i][j]*x0[j];
}
}
s=(d[i]-s)/c[i][i];
y[i]=s;
if(fabs(x0[i]-s)>temp)
{
temp=fabs(x0[i]-s);
}
}
if(temp<eps)
{
cout<<"迭代成功!迭代结果为："<<endl;
for(i=0;i<a;i++)
{
cout<<"y["<<i<<"] ="<<y[i]<<endl;
}
break;
}
if(k==m)
{
cout<<"迭代失败!!"<<endl;
break;
}
k+=1;
for(i=0;i<a;i++)
{
x0[i]=y[i];
}
}
}
void Gauss(double **c,double *d,int n,double eps)
{
int k,i;
double *y=new double[a],*x=new double[a],s,temp=0.0;
for(i=0;i<a;i++)
{
x[i]=x0[i];
y[i]=x[i];
}
k=1;
while(1)
{
temp=0.0;
for(i=0;i<a;i++)
{
s=0.0;
for(int j=0;j<a;j++)
{
if(j!=i)
{
s+=c[i][j]*y[j];
}
}
s=(d[i]-s)/c[i][i];
y[i]=s;
if(fabs(x[i]-s)>temp)
{
temp=fabs(x[i]-s);
}
}
if(temp<eps)
{
cout<<"迭代成功!迭代结果为："<<endl;
for(i=0;i<a;i++)
{
cout<<"y["<<i<<"] ="<<y[i]<<endl;
}
break;
}
if(k==m)
{
cout<<"迭代失败!!"<<endl;
break;
}
k+=1;
for(i=0;i<a;i++)
{
x[i]=y[i];
}
}
}