Task
即已知 y 0 = 0 , y 100 = 1 y_0=0,y_{100}=1 y0=0,y100=1，解线性方程组 A y = b \mathbf{A}\mathbf{y} = \mathbf{b} Ay=b，其中
A 99 × 99 = [ − ( 2 ϵ + h ) ϵ + h 0 ⋯ 0 ϵ − ( 2 ϵ + h ) ϵ + h ⋯ 0 0 ϵ − ( 2 ϵ + h ) ⋯ 0 ⋮ ⋮ ⋱ ⋱ ⋮ 0 0 ⋯ ϵ − ( 2 ϵ + h ) ] \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} A99×99= −(2ϵ+h)ϵ0⋮0ϵ+h−(2ϵ+h)ϵ⋮00ϵ+h−(2ϵ+h)⋱⋯⋯⋯⋯⋱ϵ000⋮−(2ϵ+h)
y = [ y 1 y 2 y 3 ⋮ y 99 ] b = [ a h 2 a h 2 ⋮ a h 2 a h 2 − ϵ − h ] \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*} y= y1y2y3⋮y99 b= ah2ah2⋮ah2ah2−ϵ−h
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迭代法
X ( k + 1 ) = − ( D + L ) − 1 U X ( k ) + ( D + L ) − 1 b \mathbf{X^{(k+1)}}=-(\mathbf{D}+\mathbf{L})^{-1}\mathbf{U}\mathbf{X^{(k)}}+(\mathbf{D}+\mathbf{L})^{-1}\mathbf{b} X(k+1)=−(D+L)−1UX(k)+(D+L)−1b
令
S = − ( D + L ) − 1 U = [ 0 ϵ + h 2 ϵ + h 0 ⋯ 0 0 ϵ ( ϵ + h ) ( 2 ϵ + h ) 2 ϵ + h 2 ϵ + h ⋯ 0 0 ϵ 2 ( ϵ + h ) ( 2 ϵ + h ) 3 ϵ ( ϵ + h ) ( 2 ϵ + h ) 2 ⋱ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 ϵ n − 1 ( ϵ + h ) ( 2 ϵ + h ) n 0 ⋯ ϵ ( ϵ + h ) ( 2 ϵ + h ) 2 ] \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*} S=−(D+L)−1U= 000⋮02ϵ+hϵ+h(2ϵ+h)2ϵ(ϵ+h)(2ϵ+h)3ϵ2(ϵ+h)⋮(2ϵ+h)nϵn−1(ϵ+h)02ϵ+hϵ+h(2ϵ+h)2ϵ(ϵ+h)⋮0⋯⋯⋱⋱⋯000⋮(2ϵ+h)2ϵ(ϵ+h)
I n v = ( D + L ) − 1 = [ − 1 2 ϵ + h 0 0 ⋯ 0 − ϵ ( 2 ϵ + h ) 2 − 1 2 ϵ + h 0 ⋯ 0 − ϵ 2 ( 2 ϵ + h ) 3 − ϵ ( 2 ϵ + h ) 2 − 1 2 ϵ + h ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ − ϵ n − 1 ( 2 ϵ + h ) n − ϵ n − 2 ( 2 ϵ + h ) n − 1 ⋯ − ϵ ( 2 ϵ + h ) 2 − 1 2 ϵ + h ] \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*} Inv=(D+L)−1= −2ϵ+h1−(2ϵ+h)2ϵ−(2ϵ+h)3ϵ2⋮−(2ϵ+h)nϵn−10−2ϵ+h1−(2ϵ+h)2ϵ⋮−(2ϵ+h)n−1ϵn−200−2ϵ+h1⋮⋯⋯⋯⋯⋱−(2ϵ+h)2ϵ000⋮−2ϵ+h1
在 ∥ X ( k + 1 ) − X ( k ) ∥ ∞ ≤ 1 0 − 6 \|\mathbf{X^{(k+1)}}-\mathbf{X^{(k)}}\|_{\infty}\leq10^{-6} ∥X(k+1)−X(k)∥∞≤10−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);
}
