共轭梯度法c语言实现
时间: 2024-02-05 15:03:53 浏览: 26
共轭梯度法的C语言实现如下所示:
```c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 3
void conjugateGradient(double A[N][N], double b[N], double x[N]) {
double r[N], p[N], Ap[N];
double alpha, beta, residual, residual_new;
// 初始化解向量x和残差向量r
for (int i = 0; i < N; i++) {
x[i] = 0;
r[i] = b[i];
p[i] = r[i];
}
// 计算残差向量的负梯度方向p
for (int k = 0; k < N; k++) {
// 计算矩阵A与向量p的乘积
for (int i = 0; i < N; i++) {
Ap[i] = 0;
for (int j = 0; j < N; j++) {
Ap[i] += A[i][j] * p[j];
}
}
// 计算步长alpha
double pAp = 0;
for (int i = 0; i < N; i++) {
pAp += p[i] * Ap[i];
}
alpha = pAp / (p[0] * Ap[0] + p[1] * Ap[1] + p[2] * Ap[2]);
// 更新解向量x
for (int i = 0; i < N; i++) {
x[i] += alpha * p[i];
}
// 计算新的残差向量r
for (int i = 0; i < N; i++) {
r[i] -= alpha * Ap[i];
}
// 计算残差的模
residual = 0;
for (int i = 0; i < N; i++) {
residual += r[i] * r[i];
}
residual = sqrt(residual);
// 判断是否满足收敛条件
if (residual < 1e-6) {
break;
}
// 计算beta
double r_new[N];
for (int i = 0; i < N; i++) {
r_new[i] = r[i];
}
for (int i = 0; i < N; i++) {
r_new[i] -= alpha * Ap[i];
}
double beta = 0;
for (int i = 0; i < N; i++) {
beta += r_new[i] * r_new[i];
}
beta /= residual * residual;
// 更新负梯度方向p
for (int i = 0; i < N; i++) {
p[i] = r_new[i] + beta * p[i];
}
// 更新残差
residual_new = 0;
for (int i = 0; i < N; i++) {
residual_new += r_new[i] * r_new[i];
}
residual_new = sqrt(residual_new);
// 判断是否满足收敛条件
if (residual_new < 1e-6) {
break;
}
}
}
int main() {
double A[N][N] = {{4, -1, 1}, {-1, 4, -2}, {1, -2, 4}};
double b[N] = {5, -6, 7};
double x[N];
conjugateGradient(A, b, x);
printf("Solution:\n");
for (int i = 0; i < N; i++) {
printf("%f\n", x[i]);
}
return 0;
}
```
相关推荐
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)