数值计算之二：迭代法求线性方程的解

迭代法求线性方程的解
1. 实验目标

1.设计简单迭代、高斯、SOR算法

2.对比不同方法的计算效率

2. 程序设计

本实验用C#编写，简单实现了雅可比、高斯-赛德尔迭代法、SOR迭代法等方法三种方法。三种方法抽象为ISolution接口，即有迭代次数ItNum、迭代终止条件Diff和解法Solve。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace IterEquation
{
///<summary>
///解法接口
///</summary>
interface
ISolution
{
///<summary>
///误差
///</summary>
List<double>Error {get; }
///<summary>
///迭代次数
///</summary>
int ItNum {
get; set; }

///<summary>
///终止误差
///</summary>
double Diff {
get; set; }
///<summary>
///解法
///</summary>
///<param name="A"></param>
///<param name="B"></param>
///<returns></returns>
double[] Solve(double[,]A,double[] B,double[]x0);
}
}
雅可比迭代法
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace IterEquation
{
///<summary>
///雅可比迭代方法
///</summary>
class
Jacobi:ISolution
{
private
List<double> _error =
newList<double>();
public
List<double> Error
{
get
{
return _error;
}
}

private
int _itNum =10;
public
int ItNum
{
get
{
return _itNum;
}
set
{
_itNum = value;
}
}
private
double _diff= 0.0001;
public
double Diff
{
get
{
return _diff;
}

set
{
_diff = value;
}
}

public
double[]Solve(double[,] A,
double[]B,double[] x0)
{
int it = 0;
while (it < _itNum)
{
double[] x =
newdouble[B.Length];
for (inti = 0; i < x.Length; i++)
{
doublesum = 0;

for(int j = 0; j < x.Length; j++)
{
if (i != j)
{
sum += A[i, j] *x0[j];
}
}
x[i] = (B[i] - sum) / A[i,i];
}
double diff =
Util.Diff(x0,x);

_error.Add(diff);
//误差判断
if (diff < _diff)
break;

x0 = x;
it++;
}

return x0;
}


}
}
2.高斯赛德尔迭代法

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace IterEquation
{
///<summary>
///高斯塞德尔求解法
/// GaussSeidel
///</summary>
class
GaussSeidel:ISolution
{
private
List<double> _error =
newList<double>();
public
List<double> Error
{
get
{
return _error;
}
}

private
int _itNum =10;
public
int ItNum
{
get
{
return _itNum;
}
set
{
_itNum = value;
}
}
private
double _diff= 0.0001;
public
double Diff
{
get
{
return _diff;
}

set
{
_diff = value;
}
}

public
double[]Solve(double[,] A,
double[]B, double[] x0)
{
int it = 0;
while (it < _itNum)
{
double diff = 0;
for (inti = 0; i < x0.Length; i++)
{
doublesum = 0;

for(int j = 0; j < x0.Length; j++)
{
if(i != j)
{
sum += A[i, j] *x0[j];
}
}

doublenewVal = (B[i] - sum) / A[i, i];
diff += Math.Sqrt((newVal - x0[i])*(newVal - x0[i]));
x0[i] = (B[i] - sum) / A[i,i];
}

_error.Add(diff);
if (diff < _diff)
break;

it++;
}

return x0;
}
}
}
3. SOR法

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace IterEquation
{
///<summary>
///松弛法
///</summary>
class
SOR:ISolution
{
private
List<double> _error =
newList<double>();
public
List<double> Error
{
get
{
return _error;
}
}

private
int _itNum =10;
public
int ItNum
{
get
{
return _itNum;
}
set
{
_itNum = value;
}
}
private
double _diff= 0.0001;
public
double Diff
{
get
{
return _diff;
}

set
{
_diff = value;
}
}

///<summary>
///松弛系数
///</summary>
private
double _omega= 1;

public
double Omega
{
get {
return_omega; }
set { _omega =
value;}
}

public
double[]Solve(double[,] A,
double[]B, double[] x0)
{
int it = 0;
while (it < _itNum)
{
double diff = 0;
for (inti = 0; i < x0.Length; i++)
{
doublesum = 0;

for(int j = 0; j < x0.Length; j++)
{
if (i != j)
{
sum += A[i, j] *x0[j];
}
}

doublenewVal = x0[i]+_omega*((B[i] - sum) / A[i, i] - x0[i]);
diff += Math.Sqrt((newVal - x0[i]) * (newVal - x0[i]));
x0[i] = newVal
;
}

_error.Add(diff);
if (diff < _diff)
break;

it++;
}

return x0;
}
}
}
4. 窗体程序

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using System.Windows.Forms.DataVisualization.Charting;
using System.IO;

namespace IterEquation
{
public
partial class
SolutionForm: Form
{
double[,] A;
double[] B;
double[] x0;
public SolutionForm()
{
InitializeComponent();
}

private
voidJacobiButton_Click(object sender,
EventArgs e)
{
Jacobi jsolu =
new Jacobi();
jsolu.ItNum = Int32.Parse(itnum.Text);
jsolu.Diff = Double.Parse(error.Text);

double[] x = jsolu.Solve(A, B, x0);
axisDraw(jsolu.Error, "雅可比");
resUpdate(x);
}

private
voidGaussSeidel_Click(object sender,
EventArgs e)
{
GaussSeidel gsolu =new
GaussSeidel();
gsolu.ItNum = Int32.Parse(itnum.Text);
gsolu.Diff = Double.Parse(error.Text);

double[] x = gsolu.Solve(A, B, x0);
axisDraw(gsolu.Error, "高斯-塞德尔");
resUpdate(x);

}

private
voidSOR_Click(object sender,
EventArgs e)
{
SOR gsolu =
newSOR();
gsolu.ItNum = Int32.Parse(itnum.Text);
gsolu.Diff = Double.Parse(error.Text);
gsolu.Omega = Double.Parse(omega.Text);

double[] x = gsolu.Solve(A, B, x0);
axisDraw(gsolu.Error, "SOR");
resUpdate(x);

}

///<summary>
///绘制坐标系
///</summary>
///<param name="x"></param>
///<param name="name"></param>
private
voidaxisDraw(List<double>x,String name)
{
String newName = name;

int k = 1;
while (axis.Series.FindByName(newName)!=null)
{
newName = name + k;
k++;
}

axis.Series.Add(new
Series(newName));


for(inti = 0;i < x.Count;i++)
{
axis.Series[newName].Points.AddXY(i, x[i]);
}

axis.Series[newName].ChartType =
SeriesChartType.Line;
}


///<summary>
///结果更新
///</summary>
private
voidresUpdate(double[] res)
{
String lableS =
"";
for(inti = 0;i < res.Length;i++)
{
lableS += "x" + i +"=" + res[i] +
"\n";
}
resLable.Text = lableS;
}

private
voidClearInput_Click(object sender,
EventArgs e)
{
axis.Series.Clear();
}

private
voidequation_Click(object sender,
EventArgs e)
{
ReadData();
}

private
voidReadData()
{
try
{
StreamReader objReader =new
StreamReader("equationData.txt");
string sLine =
"";
objReader.ReadLine();
sLine = objReader.ReadLine();
int numOfX =
Int32.Parse(sLine);

A = new
double[numOfX,numOfX];
B = new
double[numOfX];
x0 = new
double[numOfX];
//输入A
String[] s;
for (inti = 0; i < numOfX; i++)
{
sLine =objReader.ReadLine();
s = sLine.Split(',');
for(int j = 0; j < numOfX; j++)
{
A[i, j] = Double.Parse(s[j]);
}
}
//输入B
objReader.ReadLine();
sLine = objReader.ReadLine();
s = sLine.Split(',');
for (inti = 0; i < numOfX; i++)
{
B[i] = Double.Parse(s[i]);
}

objReader.ReadLine();
sLine = objReader.ReadLine();
s = sLine.Split(',');
for (inti = 0; i < numOfX; i++)
{
x0[i] = Double.Parse(s[i]);
}
objReader.Close();
//表示
String InputText ="";

for (inti = 0; i < numOfX; i++)
{
for(int j = 0; j < numOfX; j++)
{
if (A[i, j] < 0)
{
InputText += A[i,j] +
"x" + j;
}
else
{
if (i == 0)
{
InputText +=A[i, j] +
"x" + j;
}
else
{
InputText +=
"+" + A[i, j] + "x"+ j;
}
}
}
//输入等于号
InputText += "=" + B[i] +"\n";

}
//输入初值
InputText += "\n";
for (inti = 0; i < numOfX; i++)
{
InputText += "x"+ i +"=" + x0[i] +
",";
}

inputEqual.Text = InputText;
}
catch (IOExceptioner)
{
Console.WriteLine("An IO exception has been thrown!");
Console.WriteLine(er.ToString());
Console.ReadLine();
return;
}
}

private
voidSolutionForm_Load(object sender,
EventArgs e)
{
ReadData();
}
}
}
5. 实验分析

1．迭代法收敛效率对比

本文采用方程组在初值X1=X2=X3=0的条件下迭代10次进行测试，其中松弛条件1.2，其真值分别为X1=11，X2=12，X3=13其收敛情况如下图所示。



其中雅可比的解为，X1=10.999364458，X2=11.999362259，X3=12.999244634，高斯赛德尔的解为X1 = 10.999826484685, X2 = 11.9999893656629,X3 =12.99944028263， SOR的解为X1=11.000021026236，X2=12.0000737550456，X3=12.9999719110308。

2．松弛系数值设置范围

可以在此方程条件下看到高斯-赛德尔迭代法比雅可比迭代法迭代收敛更快，同时SOR法加快了高斯赛德尔迭代法的收敛速度。

针对SOR松弛系数的设置问题，从理论上来证明松弛系数应在0~2之间，部分方程必须小于等于1。在前面的条件下，分别取松弛度0.5,1,1.5,2四组数据对SOR进行求解如下图所示。



我们可以看到结果正如推论所示，在前三个值条件下是收敛的，且在1的时候（即退化为高斯-赛得尔迭代法）时收敛效果较好。当松弛系数等于2时，迭代不收敛。