1.cpp

线性方程组求根实验

#include<iostream>
#include<math.h>
using namespace std;
double a[11][11],b[11];
double l[11][11];//下对角矩阵l[11][11]
double y[11],x[11];
double X[11];//精确解
double R[11],X_x[11];//相关计算
void solve();
void error();
int main()
{
	//定义方程组
	for(int i=1;i<=10;i++)
	{
		for(int j=1;j<=10;j++)
			a[i][j]=1.0/(i+j-1);
		b[i]=1.0/i;
	}
	solve();
	error();
	//添加1e-7扰动
	for(int i=1;i<=10;i++)
		b[i]+=1e-7;
	solve();
	error();
	system("pause");
	return 0;
}

//平方根法解方程组
void solve()
{
	l[1][1]=sqrt(a[1][1]);
	for(int i=1;i<=10;i++)
		l[i][1]=a[i][1]/l[1][1];
	for(int j=2;j<=10;j++)
		for(int i=j;i<=10;i++)
		{
			double l_j_k_2=0.0,l_i_j_k=0.0;
			for(int k=1;k<=j-1;k++)
				l_j_k_2+=pow(l[j][k],2);
			l[j][j]=sqrt(a[j][j]-l_j_k_2);
			for(int k=1;k<=j-1;k++)
				l_i_j_k+=l[i][k]*l[j][k];
			l[i][j]=(a[i][j]-l_i_j_k)/l[j][j];
		}	
	y[1]=b[1]/l[1][1];
	for(int i=2;i<=10;i++)
	{
		double l_i_k=0.0;
		for(int k=1;k<=i-1;k++)
			l_i_k+=l[i][k]*y[k];
		y[i]=(b[i]-l_i_k)/l[i][i];
	}
	x[10]=y[10]/l[10][10];
	for(int i=9;i>=1;i--)
	{
		double l_k_i=0.0;
		for(int k=i+1;k<=10;k++)
			l_k_i+=l[k][i]*x[k];
		x[i]=(y[i]-l_k_i)/l[i][i];
	}
	//输出方程组的解
	cout<<"方程组的解x="<<endl;
	for(int i=1;i<=10;i++)
		cout<<x[i]<<endl;
}
//相关误差计算
void error()
{
	//解出所需向量
	for(int i=1;i<=10;i++)
	{
		double temp=0.0;
		for(int j=1;j<=10;j++)
			temp+=a[i][j]*x[j];
		R[i]=b[i]-temp;
		if(i==1)
			X[i]=1.0;
		else 
			X[i]=0.0;
		X_x[i]=X[i]-x[i];
	}
	//分别求各个向量的三种范数,分别是1-范数，2-范数，无穷范数
	double R_1=0.0,R_2=0.0,R_=fabs(R[1]);
	double X_x_1=0.0,X_x_2=0.0,X_x_=fabs(X_x[1]);
	double X_1=0.0,X_2=0.0,X_=fabs(X[1]);
	for(int i=1;i<=10;i++)
	{
		R_1+=fabs(R[i]);
		X_x_1+=fabs(X_x[i]);
		X_1+=fabs(X[i]);
		R_2+=pow(R[i],2);
		X_x_2+=pow(X_x[i],2);
		X_2+=pow(X[i],2);
		if(R_<fabs(R[i]))
			R_=fabs(R[i]);
		if(X_x_<fabs(X_x[i]))
			X_x_=fabs(X_x[i]);
		if(X_<fabs(X[i]))
			X_=fabs(X[i]);
	}
	R_2=sqrt(R_2);
	X_x_2=sqrt(X_x_2);
	X_2=sqrt(X_2);
	//输出所要的计算结果
	cout<<"1-范数："<<endl;
	cout<<"||r||="<<R_1;
	cout<<",||X-x||="<<X_x_1;
	cout<<",||X-x||/||X||="<<X_x_1/X_1<<endl;
	cout<<"2-范数："<<endl;
	cout<<"||r||="<<R_2;
	cout<<",||X-x||="<<X_x_2;
	cout<<",||X-x||/||X||="<<X_x_2/X_2<<endl;
	cout<<"无穷范数"<<endl;
	cout<<"||r||="<<R_;
	cout<<",||X-x||="<<X_x_;
	cout<<",||X-x||/||X||="<<X_x_/X_<<endl;
}