[算法]算法学习04（方程求解）

方程求解

二分法

对于任何连续函数，f(a)<0 f(b)>0 ，则在(a,b)中存在零点。
因此可以通过检查f((a+b)/2)的正负性来确定下一次迭代的区间。

推出条件即f((a+b)/2)==0，但是这是二分法逼近，难以达到浮点数的精确的值，因此只要判断区间小于无穷小即可。

const double PRECISION = 0.000000001;
typedef double (*FunctionPtr)(double);

double func(double x)
{
    return (2.0*x*x + 3.2*x - 1.8);
}
double DichotomyEquation(double a, double b, FunctionPtr f)
{
    double mid = (a + b) / 2.0;
    while((b - a) > PRECISION)
    {
        if(f(a) * f(mid) < 0.0) //应用迭代递推关系
            b = mid;
        else
        	a = mid;
        mid = (a + b) / 2.0; //更新迭代变量
    }
    return mid;
}

对于重根的方程式来说，难确定合适的区间进行二分，例如对于 f(x) = x^2-1 ，区间取(-2,2)时则会导致算法失效。

牛顿迭代法

先确定某点x0，求该点的导数，并以此为斜率过该点的斜线与x轴交与x1，则x1为比x0更逼近零点一点，再继续迭代。

因此，牛顿迭代法的递推公式为
x[n+1] = x[n] - f(x[n]) / f’(x[n])

double NewtonRaphson(FunctionPtr f, double x0) 
{
    double x = 0;
    int count = 0;
    do
    {
        double x1 = x0 - f(x0) / CalcDerivative(f, x0); //应用迭代递推关系
        if (fabs(x1 - x0) < PRECISION)
        {
            x = x1;
            break;
        }
        x0 = x1;  //更新迭代变量
        count++;
    } while (count < MAX_RUN_LOOP);
    return x;
}

求导部分

计算机求导智能通过导数定义，取无限小的数近似计算出导数结果。

double CalcDerivative(FunctionPtr f, double x)
{
    return (f(x + 0.000005) - f(x - 0.000005)) / 0.00001;
}

牛顿迭代法的收敛速度比二分逼近法要高。

线性代数方程组的求解

对于n阶方程组，若系数矩阵A非奇异矩阵，则方程组有唯一解，则可根据算法使结果收敛于唯一解。

雅可比迭代法（Jacobi）

推导过程：

可以理解为：左侧的系数矩阵化为单位矩阵，其他都移到等号右边，然后选取一组x的值，带入右侧式子，计算所得的左侧的值即为新的一组x。
注意推出条件一定是所有解的迭代差都小于预定精度时才可以退出。

// 后面都是，只了解了原理，具体代码没有仔细看过
const double PRECISION = 0.000000001;
const int VEC_N = 16;  //实际方程组的个数 n 必须小于VEC_N

void jacobi_eqn(double a[][VEC_N], double b[], int n, double x[])
{
    double xt[VEC_N];
    double max_diff = 0.0;

    do
    {
        for (int i = 0; i < n; i++)
        {
            double sum = 0.0;
            for (int k = 0; k < n; k++)
            {
                if (i != k)  //对角线元素不计算
                {
                    sum += a[i][k] * x[k];
                }
            }
            xt[i] = (b[i] - sum) / a[i][i];   //根据关系计算 xi 的新值
        }

        max_diff = calc_max_diff(xt, x, n);

        for (int j = 0; j < n; j++)
        {
            x[j] = xt[j]; //更新 x，准备下一次迭代
        }
    } while (max_diff > PRECISION);
}

高斯-塞德尔迭代法(Gauss-Seidal)

与Jacobi相比不同的是在每一行计算之后，将更新系数，使用新的系数带入并计算下一行，但是退出条件不好把控，不能一次计算出来，需要在迭代过程中有所记录。
具体公式：

const double PRECISION = 0.000000001;
const int VEC_N = 16;  //实际方程组的个数 n 必须小于 VEC_N

void gauss_seidel_eqn(double a[][VEC_N], double b[], int n, double x[])
{
    double max_diff = 0.0;

    do
    {
        max_diff = 0.0;
        for (int i = 0; i < n; i++)
        {
            double sum = 0.0;
            for (int k = 0; k < n; k++)
            {
                if (i != k)  //对角线元素不计算
                {
                    sum += a[i][k] * x[k];
                }
            }
            double xt = (b[i] - sum) / a[i][i];   //根据关系计算 xi 的新值
            if (std::fabs(xt - x[i]) > max_diff) //max_diff 只保留差值最大的
            {
                max_diff = std::fabs(xt - x[i]);
            }
            x[i] = xt;
        }
    } while (max_diff > PRECISION);
}

迭代法计算定积分

梯形面积法

图中的梯形面积就可以近似等于定积分结果

梯形面积法优化

用梯形面积计算时，ab距离过大，导致计算误差很大，因此可以将ab间分为n段减小误差。
公式如下：

算法实现如 trapezium() 函数所示，其定积分的精度取决于区间划分个数 n。

double trapezium(std::function<double(double)> func, double a, double b, int n)
{
    double step = (b - a) / n;
    double sum = (func(a) + func(b)) / 2.0;

    for (int j = 1; j < n; j++)
    {
        double xj = a + j*step;
        sum = sum + func(xj);
    }

    sum *= step;

    return sum;
}

变步长梯形公式法

因为对区间n等分难以确定结果的精度，在某些地方可能会产生较大误差，因此使用变步长梯形公式法计算来保证精度。
还是利用二分的思路，将ab区间分成小段去计算，而迭代的终点，是分割成的两个梯形面积与没有分割时计算的梯形面积的差值小于预定精度。

double vs_trapezium(std::function<double(double)> func, double a, double b, double eps)
{
    int n = 1;   //初始分割一个大梯形
    double t1 = (b - a) * (func(a) + func(b)) / 2.0; //用梯形公式计算初始梯形的面积
    double diff = eps + 1.0;
    do
    {
        n = 2 * n;    //梯形分割加倍
        double t = trapezium(func, a, b, n);  //用复合梯形公式法计算 n 个小梯形的面积之和
        diff = std::fabs(t1 - t);      //计算两次迭代的结果差值
        t1 = t; //更新迭代变量
    } while (diff >= eps);

    return t1;
}

辛普森公式法(Simpson)

与梯形面积法类似，可以使用中间 (a + b) / 2 的位置与 a 和 b 组成三个点进行插值计算得到求解的面积。
可以理解为取梯形的（左值+右值+4倍的中点值）的均值为底来计算面积

同样的，Simpson也可以进行n等分区间优化和变步长优化，与上面同理，此处暂时先不贴代码了。