文章目录
- Jacobi迭代方法介绍
- Gauss-Seidel迭代方法介绍
- 具体代码实现
- 示例题目
- 实现效果
Jacobi迭代方法介绍
Jacobi迭代法是一种简单的迭代求解方法,适用于严格对角占优矩阵。其基本思想是利用当前迭代步的已知解来更新下一个迭代步的解。在C语言实现中,我们首先需要定义系数矩阵A、常数向量b以及迭代向量x。然后,通过循环迭代,不断更新x的值,直到满足收敛条件或达到最大迭代次数。
Jacobi迭代法的优点是简单直观,但收敛速度相对较慢。特别是对于非严格对角占优的矩阵,Jacobi迭代法可能不收敛。因此,在实际应用中,我们需要根据问题的具体情况选择合适的迭代方法。
Gauss-Seidel迭代方法介绍
Gauss-Seidel迭代法是对Jacobi迭代法的一种改进。在Gauss-Seidel迭代法中,我们使用当前迭代步已更新的解来更新下一个解。这种策略可以加速收敛速度,特别是在对角线元素较大的情况下。在C语言实现中,我们需要对系数矩阵A进行适当的分解,以便在迭代过程中方便地获取已更新的解。
与Jacobi迭代法相比,Gauss-Seidel迭代法具有更快的收敛速度。但是,Gauss-Seidel迭代法并不总是收敛的,其收敛性取决于系数矩阵A的性质。因此,在使用Gauss-Seidel迭代法时,我们需要确保系数矩阵A满足一定的条件(如严格对角占优)。
具体代码实现
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "time.h"
#define N 3 //要求解的数组大小N x N
#define ERROR 1e-6 //误差要求
#define IterMaxNum 50000 //最大迭代次数
float current_error(float* last, float* now)
{
float max_error = 0;
float temp;
for (int i = 0; i < N; i++)
{
temp = *(now + i) - *(last + i);
if (temp < 0)
temp = -temp;
if (temp > max_error)
max_error = temp;
}
return max_error;
}
//雅可比迭代函数
float* Jacobi(float* A, float* b)
{
clock_t start, end;
double cpu_time_used;
// 记录开始时间
start = clock();
int i, j;
int iter_num = 0; //目前的迭代次数
float error_now = 100; //目前迭代下的误差
//生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
float* x = (float*)malloc(N * sizeof(float));
memset(x, 0, N * sizeof(float));
float* x0 = (float*)malloc(N * sizeof(float));
memset(x0, 0, N * sizeof(float));
while(error_now > ERROR)
{
//将上一次迭代的结果保存到x0中去
memcpy(x0, x, N * sizeof(float));
for (i = 0; i < N; i++)
{
float sum = 0;
for (j = 0; j < N; j++)
{
//printf("%f\t", A[i * N + j]);
if (i != j)
sum += A[i * N + j] * (x0[j]);
}
x[i] = (b[i] - sum) / A[i * N + i];
}
error_now = current_error(x, x0);
iter_num++;
printf("当前迭代次数为:【%d】,迭代误差为:【%f】\n", iter_num, error_now);
printf("迭代结果为:%f,%f,%f\n", x[0], x[1], x[2]);
if (iter_num > IterMaxNum)
{
printf("迭代次数超过最大值要求,返回最后一次迭代结果");
break;
}
}
// 记录结束时间
end = clock();
// 计算经过的CPU时间
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
// 打印Jacobi迭代的运行时间
printf("\n【Jacobi迭代运行了 %f 秒。】\n", cpu_time_used);
free(x0);
return x;
}
float* GS(float* A, float* b)
{
int i, j;
int iter_num = 0; //目前的迭代次数
float error_now = 100; //目前迭代下的误差
//生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
float* x = (float*)malloc(N * sizeof(float));
memset(x, 0, N * sizeof(float));
float* x0 = (float*)malloc(N * sizeof(float));
memset(x0, 0, N * sizeof(float));
clock_t start, end;
double cpu_time_used;
// 记录开始时间
start = clock();
while (error_now > ERROR)
{
memcpy(x0, x, N * sizeof(float));
for (i = 0; i < N; i++)
{
float sum = 0;
for (j = 0; j < N; j++)
{
//printf("%f\t", A[i * N + j]);
if (i != j)
sum += A[i * N + j] * (x[j]);
}
x[i] = (b[i] - sum) / A[i * N + i];
}
error_now = current_error(x, x0); //获得当前迭代的误差结果
iter_num++; //迭代次数加1
printf("当前迭代次数为:【%d】,迭代误差为:【%f】\n", iter_num, error_now);
printf("迭代结果为:%f,%f,%f\n", x[0], x[1], x[2]);
if (iter_num > IterMaxNum)
{
printf("迭代次数超过最大值要求,返回最后一次迭代结果");
break;
}
}
// 记录结束时间
end = clock();
// 计算经过的CPU时间
cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
// 打印Gauss-Seidel迭代的运行时间
printf("\n【Gauss-Seidel迭代运行了 %f 秒。】\n", cpu_time_used);
free(x0);
return x;
}
int main()
{
// 例子
float A[N][N] = {
{10,-1,-2},
{-1,10,-2},
{-1, -1, 5},
}; //要求解的矩阵
float b[N] = {72, 83, 42};
float* x;
// x = Jacobi((float*)A, b);
x = GS((float*)A, b);
printf("\n【最终结果为】:%f, %f, %f\n", x[0], x[1], x[2]);
free(x); //释放x
}
示例题目
以如下方程组为例,进行迭代方法的测试
实现效果
Jacobi迭代实现效果,迭代了【17】
达到了收敛
Gauss-Seidel迭代实现效果,而Gauss-Seidel只迭代了【10】
次即可达到收敛。
注意,有些示例中,Jacobi迭代是比Gauss-Seidel收敛的快
的,但更多的时候Gauss-Seidel是收敛更快
的,并且存在有时候Jacobi收敛
,而Gauss-Seidel却发散
,要根据情况而定。