|
|
楼主 |
发表于 2004-11-11 04:13:54
|
显示全部楼层
#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;
const int M=100;
const int N=20;
int Jacobi(double (*a)[4],double b[],double eps,int n);
int GussSaider(double (*a)[4],double b[],double eps,int n);
int SuccessiveOverRelaxationMethed(double (*a)[4],double b[],double eps,int n,double w);
int main()
{
double a[4][4]={5,2,3,2,2,4,1,-2,1,-3,4,3,3,2,2,8};
double b[4]={-1,5,4,-6};
int n=4;
double eps=1.0e-5;
cout<<"精度为"<<eps<<endl;
cout<<"用雅可比迭代法求解的结果"<<endl;
Jacobi (a, b,eps, n);
cout<<"用高斯赛德尔迭代法求解的结果"<<endl;
GussSaider (a,b,eps,n);
cout<<"用SOR迭代法求解"<<endl;
double v=1.1784;
cout<<"松弛因子为"<<v<<endl;
SuccessiveOverRelaxationMethed(a,b,eps,n,v);
return 0;
}
int Jacobi(double (*a)[4],double b[],double eps,int n)
{
int i,j,u;
double nx[N]={0},x[N]={0},d[N]={0};
for(i=0;i<n;i++)
for (j=0;j<n;j++)
if(i==j)
d=a[j];
for (u=0;u<M;u++)
{
for (i=0;i<n;i++)
{
double s=0.;
for (j=0;j<n;j++)
{
if (i==j) continue;
else
s+=a[j]*x[j];
}
nx=(b-s)/d;
}
double err=0.0;
for (j=0;j<n;j++)
if(err<fabs(nx[j]-x[j]))
err=fabs(nx[j]-x[j]);
for (j=0;j<n;j++)
x[j]=nx[j];
if (err<eps)
{ cout<<"求得的解为x_i"<<endl;
for(i=0;i<n;i++)
cout<<x<<'\t';
cout<<endl;
cout<<"迭代的次数为"<<u<<endl;
cout<<endl;
return 0;
}
}
cout<<"迭代100次后不能得到结果，求解失败"<<endl;
return 1;
}
int GussSaider(double (*a)[4],double b[],double eps,int n)
{
int i,j,u;
double s1,s2;
double nx[N]={0},x[N]={0},d[N]={0};
for(i=0;i<n;i++)
for (j=0;j<n;j++)
if(i==j)
d=a[j];
for (u=0;u<M;u++)
{
for (i=0;i<n;i++)
{ s1=0.0;
for (j=0;j<i;j++)
s1+=a[j]*nx[j];
s2=0.0;
for (j=i+1;j<n;j++)
s2+=a[j]*x[j];
nx=(b-s1-s2)/d;
}
double err=0.0;
for (j=0;j<n;j++)
if(err<fabs(nx[j]-x[j]))
err=fabs(nx[j]-x[j]);
for (j=0;j<n;j++)
x[j]=nx[j];
if (err<eps)
{ cout<<"求得的解为x_i"<<endl;
for(i=0;i<n;i++)
cout<<x<<'\t';
cout<<endl;
cout<<"迭代的次数为"<<u<<endl;
cout<<endl;
return 0;
}
}
cout<<"迭代100次后不能得到结果，求解失败"<<endl;
return 1;
}
int SuccessiveOverRelaxationMethed(double (*a)[4],double b[],double eps,int n,double w)
{
int i,j,u;
double s1,s2;
double nx[N]={0},x[N]={0},d[N]={0};
for(i=0;i<n;i++)
for (j=0;j<n;j++)
if(i==j)
d=a[j];
for (u=0;u<16;u++)
{
for (i=0;i<n;i++)
{ s1=0.0;
for (j=0;j<i;j++)
s1+=a[j]*nx[j];
s2=0.0;
for (j=i;j<n;j++)
s2+=a[j]*x[j];
nx=x+w*(b-s1-s2)/d;
}
double err=0.0;
for (j=0;j<n;j++)
if(err<fabs(nx[j]-x[j]))
err=fabs(nx[j]-x[j]);
for (j=0;j<n;j++)
x[j]=nx[j];
if (err<eps)
{ cout<<"求得的解为x_i"<<endl;
for(i=0;i<n;i++)
cout<<x<<'\t';
cout<<endl;
cout<<"迭代的次数为"<<u<<endl;
cout<<endl;
return 0;
}
}
cout<<"迭代100次后不能得到结果，求解失败"<<endl;
return 1;
} |
|
发表于 2004-11-11 03:13:18
center" align=center><FONT face="Times New Roman"><B>Jacobi迭代法</B><p></p></FONT></P><