静态迭代法

静态迭代法，是形如

\[\vec x^{(k+1)}=G\vec x^{(k)}+\vec c \]

的迭代法，而\(G\)，\(\vec c\)不随迭代步数而变化，所以是静态迭代法。常见的迭代法：雅可比法（Jacobi method）、高斯-塞德尔法（Gauss-Seidel method）和松弛法（successive relaxation method）。

静态迭代法收敛的充分条件是，存在算子范数\(||\cdot||\)，使得\(||G||<1\)。更容易计算的充分条件是\(G\)的最大特征值的模小于\(1\)。

雅可比法

线性方程组的第\(i\)个方程为

\[\sum_{j=1}^{n}a_{ij}x_j=b_i \]

假设除了\(x_i\)，其他的未知数都已经解开，那么会有

\[x_i=\frac{1}{a_{ii}}\left(b_i-\sum_{j=1,j\not=i}^{n}a_{ij}x_j\right) \]

迭代公式即为

\[x_i^{(k+1)}=\frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1}a_{ij}x_j^{(k)} - \sum_{j=i+1}^{n}a_{ij}x_j^{(k)}\right) \]

如果将矩阵\(A\)分解为三部分\(A=L+D+U\)，其中\(L\)为不包括对角线的下三角部分，\(D\)为对角线，\(U\)为不包括对角线的上三角部分，那么雅可比法的矩阵形式为

\[\vec x^{(k+1)}=-D^{-1}(L+U)\vec x^{(k)}+D^{-1}\vec b \]

高斯-塞德尔法

在雅可比法中，新值是用旧值代入解出的，而我们通常认为新值比旧值更接近精确解。如果\(x_1^{(k)}\)到\(x_n^{(k+1)}\)是按顺序依次计算的，那么在计算 \(x_i^{(k+1)}\)时可以用前面已经计算出来的 \(i-1\)个新值，即

\[x_i^{(k+1)}=\frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)} - \sum_{j=i+1}^{n}a_{ij}x_j^{(k)}\right) \]

高斯塞德尔法也可以写成矩阵形式

\[\vec x^{(k+1)}= -D^{-1}L\vec x^{(k+1)} -D^{-1}U\vec x^{(k)} + D^{-1}\vec b \]

或者

\[\vec x^{(k+1)}= -(D+L)^{-1}U\vec x^{(k)} + (D+L)^{-1}\vec b \]

这是采用从上到下的方式依次计算，其实也可以从下向上，就是将\(L\)与\(U\)互换一下就好。

松弛法

松弛法，是由高斯-塞德尔的方法经过线性加速得到的。松弛法的基础是逐次减少每一个未知的剩余。其实松弛法可以看作是高斯-塞德尔方法的推广。

在高斯-塞德尔方法的基础上，我们添加一个增量概念\(\Delta x_i^{(k)}\)，有

\[\Delta x_i^{(k)} = x_i^{(k+1)} - x_i^{(k)} \]

我们将高斯-塞德尔一般公式，可以得到

\[\begin{align} \Delta x_i^{(k)} &= x_i^{(k+1)} - x_i^{(k)} \\ &= \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)} - \sum_{j=i+1}^{n}a_{ij}x_j^{(k)}\right) - x_i^{(k)} \end{align} \]

对于迭代公式，变形为

\[\begin{align} x_i^{(k+1)} &= x_i^{(k)} + \Delta x_i^{(k)} \\ &= x_i^{(k)} + \left[ \frac{1}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)} - \sum_{j=i+1}^{n}a_{ij}x_j^{(k)}\right) - x_i^{(k)} \right] \end{align} \]

为了加快收敛速度，我们在增量\(\Delta x_i^{(k)}\)前面添加一个松弛因子\(\omega\)，最终公式为：

\[x_i^{(k+1)} = (1-\omega)x_i^{(k)} + \frac{\omega}{a_{ii}}\left( b_i - \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)} - \sum_{j=i+1}^{n}a_{ij}x_j^{(k)}\right) \]

从公式中可以看到，如果\(\omega=1\)，其实就是高斯-塞德尔公式。公式收敛条件为\(0<\omega<2\)，至于证明，没看懂。在\(0<\omega <1\)时，称为低松弛法（under-relaxation method）；当\(1<\omega <2\)时，称为超松弛法（over-relaxation method），这个条件下通常收敛较快。

逐次超松弛法（successive over-relaxation method，SOR）的矩阵形式如下：

\[\vec x^{(k+1)} = (D+\omega L)^{-1}((1-\omega)D-\omega U)\vec x^{(k)} + \omega(D+\omega L)^{-1}\vec b \]

参考资料

可能就是直接摘抄的参考资料1，还有删减……

代码实现与测试

这部分貌似是我自己写的= =，应该不是从哪里抄来的（直接根据上面公式写的）
代码实现，直接用标准库写的，那会儿我也挺闲的啊☆〜（ゝ。∂）
是不是就冲这点原创代码，我可以投稿到随笔(´▽｀)

#pragma once

#include <array>
#include <type_traits>
#include <iostream>

template<typename T, size_t M, size_t N> using Matrix = std::array<std::array<T, N>, M>;

template< typename T, size_t N >
void print(const std::array<T, N>& v) {
    for (int i=0; i<N; ++i)
        std::cout << v[i] << '\t';
    std::cout << std::endl;
}

template< typename T, size_t M, size_t N >
void print(const Matrix<T, M, N>& A) {
    for (int i=0; i<M; ++i)
        print(A[i]);
}

template< typename T, size_t M, size_t N >
std::array<T, N> multiply(const Matrix<T, M, N>& A, const std::array<T, N>& x) {
    std::array<T, M> b;
    b.fill(0);
    for (int i=0; i<M; ++i) {
        for (int j=0; j<N; ++j) {
            b[i] += A[i][j] * x[j];
        }
    }
    return b;
}

// 高斯-塞德尔法
template< typename T, size_t N >
std::array<T, N> solve_gs(
    const std::array<std::array<T, N>, N>& A,
    const std::array<T, N>& b,
    size_t max_iteration = 100) {
    static_assert (std::is_floating_point_v<T>, "floating type only");
    
    std::array<T, N> x;
    x.fill(0);
    
    std::array<T, N> D;
    for (int i=0; i<N; ++i)
        D[i] = 1.0 / A[i][i];
    
    T tolerance = 0;
    int iter = 0;
    for (; iter < max_iteration; ++iter) {
        tolerance = 0;
        for (int i=0; i<N; ++i) {
            T r = b[i];
            for (int j=0; j<N; ++j) {
                if (i!=j)
                    r -= A[i][j]*x[j];
            }
            T c = D[i] * r;
            T d = x[i] - c;
            x[i] = c;
            tolerance += d*d;
        }

        if (tolerance < 1e-12)
            break;
    }
    
    std::cout << "iteration = " << iter << ", tolerance = " << tolerance << std::endl;
    
    return x;
}

// 雅可比方法
template< typename T, size_t N >
std::array<T, N> solve_jacobi(
    const std::array<std::array<T, N>, N>& A,
    const std::array<T, N>& b, 
    size_t max_iteration = 100) {
    static_assert (std::is_floating_point_v<T>, "floating type only");
    
    std::array<T, N> x;
    x.fill(0);
    
    std::array<T, N> D;
    for (int i=0; i<N; ++i)
        D[i] = 1.0 / A[i][i];
    
    T tolerance = 0;
    int iter = 0;
    int max_iterator_count = 50;
    for (; iter < max_iteration; ++iter) {
        tolerance = 0;
        
        const std::array<T, N> y = x;
        for (int i=0; i<N; ++i) {
            T r = b[i];
            for (int j=0; j<N; ++j) {
                if (i!=j)
                    r -= A[i][j]*y[j];
            }
            x[i] = D[i] * r;
            tolerance += (y[i]-x[i]) * (y[i]-x[i]);
        }
        
        if (tolerance < 1e-12)
            break;
    }
    
    std::cout << "iteration = " << iter << ", tolerance = " << tolerance << std::endl;
    
    return x;
}

// 松弛法
template< typename T, size_t N >
std::array<T, N> solve_relax(
    const std::array<std::array<T, N>, N>& A, 
    const std::array<T, N>& b, 
    T omega, 
    size_t max_iteration = 100) {
    
    static_assert (std::is_floating_point_v<T>, "floating type only");
    
    std::array<T, N> x;
    x.fill(0);
    
    std::array<T, N> D;
    for (int i=0; i<N; ++i)
        D[i] = omega / A[i][i];
    
    T one_minus_omega = 1 - omega;
    T tolerance = 0;
    int iter = 0;
    int max_iterator_count = 50;
    for (; iter < max_iteration; ++iter) {
        tolerance = 0;
        for (int i=0; i<N; ++i) {
            T r = b[i];
            for (int j=0; j<N; ++j) {
                if (i!=j)
                    r -= A[i][j]*x[j];
            }
            T c = one_minus_omega * x[i] + D[i] * r;
            T d = x[i] - c;
            x[i] = c;
            tolerance += d*d;
        }

        if (tolerance < 1e-12)
            break;
    }
    
    std::cout << "iteration = " << iter << ", tolerance = " << tolerance << std::endl;
    
    return x;
}

using Vec3f = std::array<float, 3>;
using Vec4f = std::array<float, 4>;
using Vec3d = std::array<double, 3>;
using Vec4d = std::array<double, 4>;

using Matrix3f = Matrix<float, 3, 3>;
using Matrix4f = Matrix<float, 4, 4>;
using Matrix3d = Matrix<double, 3, 3>;
using Matrix4d = Matrix<double, 4, 4>;

using scalar = float;
using Vec3 = std::array<scalar, 3>;
using Vec4 = std::array<scalar, 4>;
using Matrix3 = Matrix<scalar, 3, 3>;
using Matrix4 = Matrix<scalar, 4, 4>;

下面是测试代码，上面的文件保存为iterative.hpp

#include <iostream>
#include "iterative.hpp"

using namespace std;

int main(int argc, char *argv[]) {

    constexpr Matrix3 A = {
        8.0f, -3.0f,  2.0f,
        4.0f, 11.0f, -1.0f, 
        6.0f,  3.0f, 12.0f
    };
    const Vec3 b = {20.0f, 33.0f, 36.0f};
    
    {
        cout << "Gauss-Seidel ----------------------------------------" << endl;
        Vec3 x = solve_gs(A, b);
        print(x);
    
        Vec3 b2 = multiply(A, x);
        print(b2);
    }
    
    {
        cout << "Jacobi ----------------------------------------------" << endl;
        Vec3 x = solve_jacobi(A, b);
        print(x);
        
        Vec3 b2 = multiply(A, x);
        print(b2);
    }
    
    cout << "*********** relax ***********" << endl;
    for (int i=0; i<20; ++i)
    {
        cout << "omega = " << 0.1f * i << "---------------------------" << endl;
        Vec3 x = solve_relax(A, b, 0.1f * i);
        print(x);
        
        Vec3 b2 = multiply(A, x);
        print(b2);
    }
    
    return 0;
}

posted @ 2021-03-07 10:55 薛定谔の三味阅读(316) 评论(0) 收藏举报

刷新页面返回顶部

薛定谔の匣

我是死还是活？不看你永远也不知道！

静态迭代法

静态迭代法

雅可比法

高斯-塞德尔法

松弛法

参考资料

代码实现与测试

公告