计算方法实验2：列主元消元法和Gauss-Seidel迭代法解线性方程组

Task

在这里插入图片描述
即已知 $y_0=0,y_{100}=1$ ，解线性方程组 $\mathbf{A}\mathbf{y} = \mathbf{b}$ ，其中

$\mathbf{A_{99\times99}}= \begin{bmatrix} -(2\epsilon + h) & \epsilon + h & 0 & \cdots & 0 \\ \epsilon & -(2\epsilon + h) & \epsilon + h & \cdots & 0 \\ 0 & \epsilon & -(2\epsilon + h) & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots& \epsilon& -(2\epsilon + h) \\ \end{bmatrix}$

$\begin{align*} \mathbf{y} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{99} \\ \end{bmatrix} \quad & \quad \mathbf{b}= \begin{bmatrix} ah^2 \\ ah^2 \\ \vdots \\ ah^2 \\ ah^2-\epsilon-h \end{bmatrix} \end{align*}$
在这里插入图片描述

Algorithm1：列主元消元法

Code

#include<iostream>
#include<cmath>
#include <iomanip>
using namespace std;

long double a = 1.0 / 2, n = 100, h = 1.0 / n;//注意a=1/2会把a赋为0而不是0.5
long double x[99], A[99][99], b[99];

// 交换两个数组的元素
void swap(long double* x, long double* y)
{
    for(int i = 0; i <= n - 1; i++)
    {
        long double temp = x[i];
        x[i] = y[i];
        y[i] = temp;
    }
    return ;
}

long double* Column_Elimination(long double (*A)[99])
{
    long double* x = new long double[99]();
    long double (*matrix)[100] = new long double[99][100];//增广矩阵
    fill_n(&matrix[0][0], 99*100, 0);
    for(int i = 0; i < n - 1; i++)
        for(int j = 0; j < n - 1; j++)
            matrix[i][j] = A[i][j];
    for(int i = 0; i < n - 1; i++)
        matrix[i][99] = b[i];
    for(int col = 0; col < n - 1; col++)//找到列主元
    {
        long double maxnum = abs(matrix[col][col]);
        int maxrow = col;
        for(int row = col + 1; row < n - 1; row++)
        {
            if(abs(matrix[row][col]) > maxnum)
            {
                maxnum = abs(matrix[row][col]);
                maxrow = row;
            }
        }
        swap(matrix[col], matrix[maxrow]);
        for(int row = col + 1; row < n - 1; row++)
        {
            long double res = matrix[row][col] / matrix[col][col];
            for(int loc = col; loc <= n - 1; loc++)
                matrix[row][loc] -= matrix[col][loc] * res; 
        }
    }
    for(int row = n - 2; row >= 0; row--)//回代
    {
        for(int col = row + 1; col < n - 1; col++)
        {
            matrix[row][99] -= matrix[col][99] * matrix[row][col] / matrix[col][col];
            matrix[row][col] = 0;
        }
        matrix[row][99] /= matrix[row][row];
        matrix[row][row] = 1;
    }
    for(int i = 0; i < n - 1; i++)
        x[i] = matrix[i][99];
    return x;
}

long double* PreciseSol(long double* x , long double epsilon)
{
    long double* y = new long double[99];
    for(int i = 0; i < n - 1; i++)
        y[i] = (1 - a) * (1 - exp(-x[i] / epsilon)) / (1 - exp(-1 / epsilon)) + a * x[i];
    return y;
}

void Calculate(long double epsilon)
{
    for(int i = 0; i < n - 1; i++)
        b[i] = a * h * h;
    b[98] -= epsilon + h;
    long double* Presol = PreciseSol(x , epsilon);
    for(int i = 0; i < n - 1; i++)
    {
        A[i][i + 1] = epsilon + h;
        A[i + 1][i] = epsilon;
    }
    for(int i = 0; i < n - 1; i++)
        A[i][i] = -(2 * epsilon + h);
    long double* Column = Column_Elimination(A);
    long double errorColumn = 0;
    for(int i = 0; i < n - 1; i++){
        errorColumn += abs(Column[i] - Presol[i]) / 99;
        //cout << Column[i] << " " << Seidel[i] << " " << Presol[i] << endl;
    }
    cout << "epsilon=" << epsilon << "时，列主元消元法的相对误差为"  << errorColumn << endl;
    delete[] Presol;
}

int main()
{
    for(int i = 0; i < n; i++)
        x[i] = (i + 1) * h;
    Calculate(1);
    Calculate(0.1);
    Calculate(0.01);
    Calculate(0.0001);
}

Algorithm2：Gauss-Seidel迭代法

$\mathbf{X^{(k+1)}}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}\mathbf{X^{(k)}}+(\mathbf{D}+\mathbf{L})^{-1}\mathbf{b}$
令
$\begin{align*} \mathbf{S} &= -(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U} \\ &= \begin{bmatrix} 0 & \frac{\epsilon + h}{2\epsilon + h} & 0 & \cdots & 0 \\ 0 & \frac{\epsilon(\epsilon + h)}{(2\epsilon + h)^2} & \frac{\epsilon + h}{2\epsilon + h}& \cdots & 0 \\ 0 & \frac{\epsilon^2(\epsilon + h)}{(2\epsilon + h)^3} & \frac{\epsilon(\epsilon + h)}{(2\epsilon + h)^2} & \ddots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \frac{\epsilon^{n-1}(\epsilon + h)}{(2\epsilon + h)^n} & 0 & \cdots & \frac{\epsilon(\epsilon + h)}{(2\epsilon + h)^2} \\ \end{bmatrix} \end{align*}$
$\begin{align*} \mathbf{Inv} &= (\mathbf{D}+\mathbf{L})^{-1} \\ &= \begin{bmatrix} -\frac{1}{2\epsilon + h} & 0 & 0 & \cdots & 0 \\ -\frac{\epsilon}{(2\epsilon + h)^2} & -\frac{1}{2\epsilon + h} & 0 & \cdots & 0 \\ -\frac{\epsilon^2}{(2\epsilon + h)^3} & -\frac{\epsilon}{(2\epsilon + h)^2} & -\frac{1}{2\epsilon + h} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ -\frac{\epsilon^{n-1}}{(2\epsilon + h)^n} & -\frac{\epsilon^{n-2}}{(2\epsilon + h)^{n-1}} & \cdots & -\frac{\epsilon}{(2\epsilon + h)^2} & -\frac{1}{2\epsilon + h} \\ \end{bmatrix} \end{align*}$
在 $\|\mathbf{X^{(k+1)}}-\mathbf{X^{(k)}}\|_{\infty}\leq10^{-6}$ 时结束迭代。

Code

#include<iostream>
#include<cmath>
#include <iomanip>
using namespace std;

long double a = 1.0 / 2, n = 100, h = 1.0 / n;
long double x[99], A[99][99], b[99];

// 矩阵和向量的乘法，用于求Gauss-Seidel迭代的f向量
long double* Multiplie(long double (*Inv)[99], long double* x)
{
    long double* res = new long double[99]();
    for(int i = 0; i < n - 1; i++)
    {
        for(int j = 0; j < n - 1; j++)
            res[i] += Inv[i][j] * x[j];
    }
    return res;
}

// 计算两个向量的差向量的无穷范数
long double Norm(long double* x1, long double* x2)
{
    long double norm = 0;
    for(int i = 0; i < n - 1; i++)
        if(abs(x1[i] - x2[i]) > norm)
            norm = abs(x1[i] - x2[i]);
    return norm;
}

// Gauss-Seidel迭代法求解线性方程组
long double* Gauss_Seidel(long double (*S)[99], long double (*Inv)[99])
{
    long double* x1 = new long double[99]();
    long double tempx[99];
    for(int i = 0; i < n - 1; i++)
        tempx[i] = x[i];
    long double *temp1 = Multiplie(S, tempx);
    long double *temp2 = Multiplie(Inv, b);
    for(int i = 0; i < n - 1; i++)    
        x1[i] = temp1[i] + temp2[i];
    while(Norm(tempx, x1) > 1e-6)
    {
        for(int i = 0; i < n - 1; i++) 
            tempx[i] = x1[i];
        long double *temp1 = Multiplie(S, tempx);
        long double *temp2 = Multiplie(Inv, b);
        for(int i = 0; i < n - 1; i++)    
            x1[i] = temp1[i] + temp2[i];
    }
    return x1;
}

// 计算精确解
long double* PreciseSol(long double* x , long double epsilon)
{
    long double* y = new long double[99];
    for(int i = 0; i < n - 1; i++)
        y[i] = (1 - a) * (1 - exp(-x[i] / epsilon)) / (1 - exp(-1 / epsilon)) + a * x[i];
    return y;
}

// 计算误差并输出
void Calculate(long double epsilon)
{
    for(int i = 0; i < n - 1; i++)
        b[i] = a * h * h;
    b[98] -= epsilon + h;
    long double* Presol = PreciseSol(x , epsilon);
    long double S[99][99], Inv[99][99];
    fill(&S[0][0], &S[0][0] + sizeof(S) / sizeof(S[0][0]), 0);
    fill(&Inv[0][0], &Inv[0][0] + sizeof(Inv) / sizeof(Inv[0][0]), 0);
    long double init1 = (epsilon + h) / (2 * epsilon + h);
    long double temp = epsilon / (2 * epsilon + h);
    for(int i = 0; i < n - 1; i++)//初始化S=-(D+L)^(-1)U
    {
        for(int j = i , k = 1; j < n -1 && k < n - 1; j++, k++)
        {
            S[j][k] = init1;
        }
        init1 *= temp; 
    }
    long double init2 = -1 / (2 * epsilon + h);
    for(int i = 0; i < n - 1; i++)//初始化(D+L)^(-1)
    {
        for(int j = i, k = 0; j < n - 1 && k < n - 1; j++ , k++)
        {
            Inv[j][k] = init2;
        }
        init2 *= temp; 
    }
    long double* Seidel = Gauss_Seidel(S, Inv);
    long double errorSeidel = 0;
    for(int i = 0; i < n - 1; i++)
        errorSeidel += abs(Seidel[i] - Presol[i]) / 99;
    cout << "epsilon=" << epsilon << "时，Gauss_Seidel迭代法的相对误差为"  << errorSeidel << endl;
    delete[] Presol;
}

int main()
{
    for(int i = 0; i < n - 1; i++)
        x[i] = (i + 1) * h;
    Calculate(1);
    Calculate(0.1);
    Calculate(0.01);
    Calculate(0.0001);
}