追赶法(Thomas) 雅克比迭代(Jacobi) 高斯迭代(Gauss) 的C++实现

River Chandler

已于 2022-12-19 21:44:58 修改

阅读量615

点赞数 1

分类专栏： C系列语言程序设计方法与案例 # 科学计算原理文章标签： C++ 算法 python

于 2022-12-18 00:00:06 首次发布

本文链接：https://blog.csdn.net/Chandler_river/article/details/128358079

版权

科学计算原理同时被 2 个专栏收录

24 篇文章 2 订阅

订阅专栏

C系列语言程序设计方法与案例

5 篇文章 0 订阅

订阅专栏

其实如果你搜到了这个博客
- 基本可以说明,你肯定已经知道什么是追赶法了
- Python Python 就是大一大二的小孩拿来使自我感觉良好的骗子玩具(对低龄儿童和专业人士除外)

$\bg_white \begin{bmatrix} b_1 &c_1 & & & \\ a_2& b_2 &c_2 & & \\ & \ddots \ddots& \ddots & & \\ & & a_{n-1} & b_{n-1} & c_{n-1}\\ & & & a_n& b_n \end{bmatrix}\begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_4\\ x_5 \end{bmatrix}=\begin{bmatrix} d_1\\ d_2\\ \vdots \\ d_{n-1}\\ d_n \end{bmatrix}$

LU分解

$\bg_white \begin{bmatrix} b_1 &c_1 & & & \\ a_2& b_2 &c_2 & & \\ & \ddots \ddots& \ddots & & \\ & & a_{n-1} & b_{n-1} & c_{n-1}\\ & & & a_n& b_n \end{bmatrix}=\begin{bmatrix} \alpha_1 & & & & \\ \gamma_2& \alpha_2 & & & \\ & \ddots & \ddots & & \\ & & \gamma_{n-1}& \alpha_{n-1} & \\ & & & \gamma_n&\alpha_n \end{bmatrix}\begin{bmatrix} 1 &\beta_1 & & & \\ & 1 & \beta_2 & & \\ & & \ddots &\ddots & \\ & & & 1 &\beta_{n-1} \\ & & & & 1 \end{bmatrix}$

$\left\{\begin{matrix} \gamma_i=a_i\,\,\,\,\,\,\,\,\,\,(i=2,3,4,...,n)\\ b_i=\gamma_i\beta_{i-1}+\alpha_i\,\,\,(i=2,3,4,...,n) \\ c_i=\alpha_i \beta_i\,\,\,\,\,\,\,\,(i=1,2,3,...,n-1) \end{matrix}\right.$

LU分解后,可以轻松得得到每个方程得解
哪两个方程呢:
- AX=b
- LUX=b====>LY=b====>Y=d
- UX=d====>X

$\begin{matrix} y_1=d_1/b_1,y_i=(d_i-a_iy_{i-1})/(b_i-a_i\beta_{i-1}),i=2,3,...,n\\ x_n=y_n,x_i=y_i-\beta_ix_{i+1},i=n-1,n-2,...,1 \end{matrix}$

C++ 实现
其实有个很有趣的事情,...,这玩意不就是顺下来的...
啊,对对对

#include <iostream>
#include <vector>
using namespace std;

vector<double> THOMAS(vector<vector<double>>a, vector<double>d)
{
	int n = d.size();
	
	vector<double>x(n);
	vector<double>alpha(n);
	vector<double>gamma(n);
	vector<double>beta(n);
	vector<double>y(n);
	
	alpha[0] = a[0][0];
	beta[0] = a[0][1] / a[0][0]; y[0] = d[0] / a[0][0];
	
	for (int i=1;i<n;i++)
	{
		gamma[i] = a[i][i-1];
		alpha[i] = a[i][i] - gamma[i] * beta[i-1];
		if (i<n-1)
		{
			beta[i] = a[i][i+1]/alpha[i];
		}
		y[i] = (d[i]-a[i][i-1]*y[i-1])/alpha[i];
	}
	x[n-1] = y[n-1];
	for (int i=n-2;i>=0;i--)
	{
		x[i] = y[i]-beta[i]*x[i+1];
	}
	return x;
}
int main()
{
	vector<vector<double>>a(4,vector<double>(4));
	vector<double>b(4);
	vector<double>x(4);
	
	a[0] = {2,1,0,0}; 
	a[1] = {1,2,1,0}; 
	a[2] = {0,1,2,1}; 
	a[3] = {0,0,1,2};
	b    = {1,1,1,1};
	
	x = THOMAS(a, b);
	
	for (int i=0;i<4;i++)
	{
		for (int j=0;j<4;j++)
		{    
			cout << a[i][j] << " ";
		}
		cout << b[i] << endl;
	}
	cout << "SOLUTION of THOMAS method:" << endl;
	for (int i=0;i<4;i++)
	{
		cout<<"x"<<i+1<<" = " <<x[i] << endl;
	}
}

显然,这是我们所不能满意的

我们希望的功能有

1.检验是否可以使用THOMAS方法

2.可以从文件中读取,好,明天更

3.明天cupt 裁判