1: #ifndef CHASEMETHORD_H
2: #define CHASEMETHORD_H
3: #include <assert.h>
4: #define NULL 0
5: /************************************************************************/
6: /* 功 能：追赶法接矩阵方程：Ax = y.
7: /* 其中A必须是三对角矩阵，即非零元素分布在主对角线及其相邻两条次对角
8: /* 线上。
9: /* 参 数：size_t n, 矩阵对角线元素个数[in]
10: /* T* a, 主对角线以下元素[in]
11: /* T* b, 主对角线以上元素[in]
12: /* T* c, 主对角线元素[in]
13: /* T* x, 矩阵的解集[out]
14: /* T* y, y集[in]
15: /************************************************************************/
16: template <typename T>
17: void ChaseMethord(size_t n, T* a, T *b, T* c, T* x, T* f)
18: { 19: assert(a != NULL && b != NULL && c != NULL && x != NULL && f != NULL);
20:
21: //Ax = f => A = LU,即L(Ux) = f,又：令y = Ux,所以Ly = f;
22: //其中，L是下二对角矩阵，U是单位上二对角矩阵
23: T *l = new T[n];
24: T *u = new T[n];
25: T *y = new T[n];
26:
27: assert(l != NULL);
28: assert(u != NULL);
29: assert(y != NULL);
30:
31: l[0] = b[0];
32: if (l[0] == 0)
33: { 34: return;
35: }
36: y[0] = f[0] / l[0];
37:
38: size_t i;
39:
40: //依次求出l[0]->u[0]->l[1]->u[1]->...->l[n-1],y[i]
41: for (i = 1; i < n; i++)
42: { 43: if (l[i] == 0)
44: { 45: return;
46: }
47:
48: u[i - 1] = c[i - 1] / l[i - 1];
49: l[i] = b[i] - a[i - 1] * u[i - 1];
50: y[i] = (f[i] - a[i - 1] * y[i - 1]) / l[i];
51: }
52:
53: //逆序求x[i]
54: x[i-1] = y[i-1];
55: for (int k = i - 2; k >= 0; k--)
56: { 57: x[k] = y[k] - u[k] * x[k + 1];
58: }
59:
60: if (l != NULL)
61: { 62: delete [] l;
63: l = NULL;
64: }
65:
66: if (u != NULL)
67: { 68: delete [] u;
69: u = NULL;
70: }
71:
72: if (y != NULL)
73: { 74: delete [] y;
75: y = NULL;
76: }
77: }
78:
79:
80: #endif
1: #include "ChaseMethord.h"
2: #include <iostream>
3: using namespace std;
4:
5: int main()
6: { 7: float a[] = {1,1,2}; 8: float b[] = {2,3,1,1}; 9: float c[] = {1,1,1}; 10: float f[] = {1,2,2,0}; 11:
12: float *x = new float[4];
13:
14: ChaseMethord(4, a, b, c, x,f);
15:
16: for (int i = 0; i < 4; i ++)
17: { 18: cout << x[i] << endl;
19: }
20:
21: delete [] x;
22:
23: return 0;
24: }