线性方程组的迭代解法
class Operation
{
double[][] a; //系数矩阵
double[] b; //常数向量
double[] x00,x0, x; //初始迭代向量和迭代向量
int n; //矩阵大小
public double epsilon = 1e-5; //最大精度
int M, m; //最大迭代向量次数和迭代次数
double current_e; //当前最大精度值
public Operation(int n)
{
this.n = n;
a = new double[n][];
for(int i = 0; i < n;i++)
{
a[i] = new double[n];
}
b = new double[n];
x0 = new double[n];
x = new double[n];
m = 0;
}
public void Init(double[][] a, double[] b,double[] x00, double epsilon = 1e-5, int M = 5000) //初始化
{
this.a = a;
this.b = b;
this.epsilon = epsilon;
this.M = M;
this.x00 = x00;
for (int i = 0; i < n; i++)
x0[i] = x00[i];
}
public void Change_x0(double[] x0) //改变初始向量
{
this.x00 = x0;
for (int i = 0; i < n; i++)
this.x0[i] = x00[i];
}
public void Change_e(double e) //改变精度
{
epsilon = e;
}
public void ReSet() //重置向量
{
for (int i = 0; i < n; i++)
x0[i] = x00[i];
x = new double[n];
m = 0;
current_e = 0;
}
public double[] Result() //返回迭代向量
{
return x;
}
public int Reslut_m() //返回迭代次数
{
return m;
}
public void Show()
{
Console.Write("迭代向量:");
for (int i = 0; i < 3; i++)
Console.Write("x" + (i + 1) + "=" + Result()[i] + " ");
Console.WriteLine();
Console.WriteLine("迭代次数:" + Reslut_m());
Console.WriteLine();
ReSet();
}
public void Jacobi() //雅克比迭代法
{
do
{
current_e = 0;
for(int i = 0;i < n;i++)
{
double sum = 0;
for (int j = 0; j < n; j++)
{
if (j != i)
{
sum = sum + a[i][j] * x0[j];
}
}
x[i] = (b[i] - sum) / a[i][i];
} //更新迭代向量
m++; //迭代次数+1
for (int i = 0; i < n; i++)
{
if (Math.Abs(x[i] - x0[i]) > current_e)
current_e = Math.Abs(x[i] - x0[i]);
} //计算当前误差
for (int i = 0; i < n; i++)
x0[i] = x[i]; //更新初始向量
} while (current_e > epsilon && m < M); //判断是否仍未达到精度要求且未达到最大迭代次数
}
public void Gauss_Seidel() //高斯-赛德尔迭代法
{
do
{
current_e = 0;
for (int i = 0; i < n; i++)
{
double sum = 0;
for (int j = 0; j < n; j++)
{
if (j != i)
{
if(j < i)
{
sum = sum + a[i][j] * x[j];
}
else
{
sum = sum + a[i][j] * x0[j];
}
}
}
x[i] = (b[i] - sum) / a[i][i];
} //更新迭代向量
m++; //迭代次数+1
for (int i = 0; i < n; i++)
{
if (Math.Abs(x[i] - x0[i]) > current_e)
current_e = Math.Abs(x[i] - x0[i]);
} //计算当前误差
for (int i = 0; i < n; i++)
x0[i] = x[i]; //更新初始向量
} while (current_e > epsilon && m < M); //判断是否仍未达到精度要求且未达到最大迭代次数
}
public void SOR(double omega) //SOR迭代算法
{
do
{
current_e = 0;
for (int i = 0; i < n; i++)
{
double sum = 0;
for (int j = 0; j < n; j++)
{
if (j != i)
{
if (j < i)
{
sum = sum + a[i][j] * x[j];
}
else
{
sum = sum + a[i][j] * x0[j];
}
}
}
x[i] = (b[i] - sum) / a[i][i];
x[i] = omega * x[i] + (1 - omega) * x0[i];
} //更新迭代向量
m++; //迭代次数+1
for (int i = 0; i < n; i++)
{
if (Math.Abs(x[i] - x0[i]) > current_e)
current_e = Math.Abs(x[i] - x0[i]);
} //计算当前误差
for (int i = 0; i < n; i++)
x0[i] = x[i]; //更新初始向量
} while (current_e > epsilon && m < M); //判断是否仍未达到精度要求且未达到最大迭代次数
}
}
用法演示
static void Main(string[] args)
{
double[][] a = { new double[] { 3, 2, 1 }, new double[] { 0, 5, 1 }, new double[] { 0, 4, 8 } };
double[] b = { 39, 24, 39 };
Operation q = new Operation(3);
q.Init(a, b, new double[] { 1, 1, 1 });
Console.WriteLine("雅克比迭代法:");
q.Jacobi();
q.Show();
q.ReSet();
Console.WriteLine("高斯赛德尔迭代法:");
q.Gauss_Seidel();
q.Show();
q.ReSet();
}