数値計算の基本をC++で書いていきます.
boostやEigenを使用しているコードのコンパイルには, boostやEigenのインストールがそれぞれ必要です.
変更した方が良い箇所があれば, この記事かTwitterの@Yonono_01まで連絡ください.
目次
連立方程式 (Simultaneous Equation)
数値計算ではそのモデルを表す方程式を離散化する.
このとき, 連立一次方程式の問題に帰着されることが多い.
連立一次方程式の解法としては, Gaussの消去法やGauss Seidel法などがある.
Gaussの消去法 (Gaussian Elimination)
Gaussの消去法は係数行列を掃き出し法により上三角行列に変形し, 後退代入により解を求めるものである.
このとき, 部分pivot選択をおこない, 丸め誤差の影響を小さくする.
#include <cmath> // for std::fabs, std::pow
#include <iostream> // for std::cerr, std::cin, std::cout
#include <optional> // for std::optional
#include <string> // for std::string
#include <vector> // for std::vector
#include <utility> // for std::move
#include <boost/format.hpp> // for boost::format
#include <boost/multi_array.hpp> // for boost::multi_array
#include <boost/range/algorithm.hpp> // for boost::range::copy
namespace{
using dmatrix = boost::multi_array<double, 2>;
using dvector = std::vector<double>;
static auto constexpr MAXBUFSIZE = 32;
static auto const eps = pow(2, -50);
template <typename T, bool check_positive_num = true>
T input_parameter(std::string const &str){
T param;
while (true){
std::cout << str;
std::cin >> param;
std::cin.clear();
std::cin.ignore(MAXBUFSIZE, '\n');
// 初期条件として<=0の場合に再入力が要らない時は<type, false>
if constexpr (check_positive_num){
if (!std::cin.fail() && param > static_cast<T>(0)){
break;
}
}
else{
if (!std::cin.fail()){
break;
}
}
}
return param;
}
void display_result(dvector const &x);
// 後退代入
void Gauss_backward(dmatrix &A);
// Gaussの消去法
std::optional<dvector> Gauss_elimination(dmatrix &A, double eps);
// 掃き出し法
bool Gauss_forward(dmatrix &A, double eps);
// pivot選択
void Gauss_pivot(dmatrix &A, int k);
dvector input_row(int size);
dvector make_result(dmatrix const &A);
}
int main(){
auto const size = input_parameter<int> ("Enter the size of a expansion coefficient matrix n-rows and (n+1)-colowns\n");
dmatrix A(boost::extents[size][size + 1]);
std::cout << "Enter elements of the expansion coefficient matrix" << std::endl;
for (int i = 0; i < size; i++){
auto const row = input_row(size + 1);
boost::range::copy(row, A[i].begin());
}
if (auto const res = Gauss_elimination(A, eps); res){
auto const x(*res);
display_result(x);
}
else{
std::cerr << "Error! Gaussian Elimination failed" << std::endl;
return -1;
}
}
namespace{
void display_result(dvector const &x){
auto const size = static_cast<int>(x.size());
for (int i = 0; i < size; i++){
std::cout << "x[" << i << boost::format("] = %.7f\n") % x[i];
}
}
void Gauss_backward(dmatrix &A){
auto const size = static_cast<int>(A.size());
A[size - 1][size] /= A[size - 1][size - 1];
for (int i = size - 2; i >= 0; i--){
auto ax = 0.0;
for (int j = i + 1; j < size; j++){
ax += A[i][j] * A[j][size];
}
A[i][size] = (A[i][size] - ax) / A[i][i];
}
}
std::optional<dvector> Gauss_elimination(dmatrix &A, double eps){
if (!Gauss_forward(A, eps)){
return std::nullopt;
}
Gauss_backward(A);
return std::make_optional(make_result(A));
}
bool Gauss_forward(dmatrix &A, double eps){
auto const size = static_cast<int>(A.size());
for (int k = 0; k < size - 1; k++) {
Gauss_pivot(A, k);
for (int i = k + 1; i < size; i++){
for (int j = k; j < size + 1; j++){
A[i][j] -= A[i][k] * (A[k][j] / A[k][k]);
auto tmp = A[i][j];
auto a = 0;
}
}
}
if (std::fabs(A[size - 1][size - 1]) < eps){
return false;
}
return true;
}
void Gauss_pivot(dmatrix &A, int k){
auto const size = static_cast<int>(A.size());
auto pivot_ren = k;
auto pivot = std::fabs(A[k][k]);
for (int i = k + 1; i < size; i++){
if (std::fabs(A[i][k]) > pivot){
pivot_ren = i;
pivot = std::fabs(A[i][k]);
}
}
if (k != pivot_ren){
dvector tmp(A[k].begin(), A[k].end());
boost::range::copy(A[pivot_ren], A[k].begin());
boost::range::copy(tmp, A[pivot_ren].begin());
}
}
dvector input_row(int size){
dvector row(size);
auto success = true;
do{
auto i = 0;
for (; i < size; i++) {
std::cin >> row[i];
if (std::cin.fail()) {
success = false;
break;
}
}
if (i == size) {
success = true;
}
std::cin.clear();
std::cin.ignore(MAXBUFSIZE, '\n');
} while (!success);
return row;
}
dvector make_result(dmatrix const &A){
auto const size = static_cast<int>(A.size());
dvector x(size);
for (int i = 0; i < size; i++){
x[i] = A[i][size];
}
return x;
}
}
Gauss Seidel法 (Gauss Seidel Method)
Gauss Seidel法は連立一次方程式の解を反復法によって求めるものである.
Gauss Seidel法は係数行列が正定値対称であるとき収束する.
また, 係数行列が狭義対角優位であるとき収束する.
#include <cassert> // for assert
#include <cmath> // for std::fabs, std::pow, std::sqrt
#include <iostream> // for std::cerr, std::cin, std::cout
#include <optional> // for std::optional
#include <string> // for std::string
#include <utility> // for std::move
#include <vector> // for std::vector
namespace{
static auto constexpr MAXBUFSIZE = 32;
static auto constexpr MAXLOOP = 100000;
static auto const eps = std::pow(2.0, -50);
template <typename T, bool check_positive_num = true>
T input_parameter(std::string const & str){
T param;
while (true){
std::cout << str;
std::cin >> param;
std::cin.clear();
std::cin.ignore(MAXBUFSIZE, '\n');
// 初期条件として<=0の場合に再入力が要らない時は<type, false>
if constexpr (check_positive_num){
if (!std::cin.fail() && param > static_cast<T>(0)){
break;
}
}
else{
if (!std::cin.fail()){
break;
}
}
}
return param;
}
double calc_norm(std::vector<double> const &afterx, std::vector<double> const &beforex);
std::optional<std::vector<double>> Gauss_Seidel(std::vector<std::vector<double>> const &A, std::vector<double> const &b, double eps);
double inner_product(std::vector<double> const &a, std::vector<double> const &b);
}
int main(){
auto const size = input_parameter<int> ("Enter the size of a coefficient matrix\n");
// Ax = b
std::vector<std::vector<double>> A(size, (std::vector<double>(size)));
std::vector<double> b(size);
std::cout << "Enter elements of the coefficient matrix" << std::endl;
for (int i = 0; i < size; i++){
for (int j = 0; j < size; j++){
A[i][j] = input_parameter<double, false>("");
}
}
std::cout << "Enter elements of a right side matrix" << std::endl;
for (int i = 0; i < size; i++){
b[i] = input_parameter<double, false>("");
}
// Gauss Seidel methodが成功して解が収束していれば, resに解が入る
if (auto const res = Gauss_Seidel(A, b, eps); res){
// resから中身を取り出す
auto const x = *res;
for (auto i = 0; i < size; i++){
std::cout << x[i] << std::endl;
}
}
else{
// resに中身が入っていなかったら、解が収束していないのでエラーメッセージ
std::cerr << "Error! Gauss_Seidel method failed!" << std::endl;
return -1;
}
}
namespace{
double calc_norm(std::vector<double> const &afterx, std::vector<double> const &beforex){
auto const size = afterx.size();
// afterxとbeforexのサイズが違っていたら強制終了
assert(size == beforex.size());
auto norm = 0.0;
for (auto i =0ULL; i < size; i++){
norm += (afterx[i] - beforex[i]) * (afterx[i] - beforex[i]);
}
return std::sqrt(norm);
}
std::optional<std::vector<double>> Gauss_Seidel(std::vector<std::vector<double>> const &A, std::vector<double> const &b, double eps){
auto const size = A.size();
// vect Aとvect bのサイズが違っていたら強制終了
assert(size == b.size());
std::vector<double> x(size, 0.0);
for (auto i = 0; i < MAXLOOP; i++){
auto const beforex(x);
for (auto j =0ULL; j < size; j++){
x[j] = (b[j] - (inner_product(A[j], x) - A[j][j] * x[j])) / A[j][j];
}
// 残差vectorのnormがEPS未満になったら終了
if (calc_norm(x, beforex) < eps){
// 収束したので、解ベクトルxをstd::moveで右辺値参照にキャストし、std::optionalでラップして返す
return std::make_optional(std::move(x));
}
}
return std::nullopt;
}
double inner_product(std::vector<double> const &a, std::vector<double> const &b){
auto value = 0.0;
auto const size = a.size();
// vect aとvect bのサイズが違っていたら強制終了
assert(size == b.size());
for (auto i = 0ULL; i < size; i++){
value += a[i] * b[i];
}
return value;
}
}
TDMA法 (Tri-Diagonal Matrix Algorithm method)
執筆中
その他の解法
連立一次方程式を解く方法としては他にはCramerの公式(Cramer rule)やJacobi法(Jacobi rule)などがある.
Gauss Seidel法における計算で緩和係数を考えて調整するものをSOR法(Successive Over-Relaxation)という.
実際の数値計算ではCG法(conjugate gradient method)やCholesky分解(Cholesky decomposition), LU分解(LU decomposition)が主に用いられる.
Eigen
C++には行列計算ライブラリとしてEigenがある.
行列に関する計算は基本的にEigenを使う.
#include <iostream> // for std::cout, std::endl;
#include <Eigen/Dense> // for Eigen::MatrixXd, Eigen::VectorXd, Eigen::ColPivHouseholderQR
int main(){
Eigen::Matrix3d A;
Eigen::Vector3d b;
A << 2,3,1, 4,-5,-1, 6,-7,1;
b << 7, 3, 5;
Eigen::ColPivHouseholderQR<Eigen::Matrix3d> dec(A);
Eigen::Vector3d x = dec.solve(b);
std::cout << x << std::endl;
}
補遺 Eigenについて
EigenはC++のテンプレートを使用されて作られている線形代数ライブラリである.
Eigenの基本的なコマンドを紹介する.
ここで紹介していない便利なコマンドもたくさんあるので, 参考文献を適宜参考にしてほしい.
#include <iostream> // for std::cout, std::endl
#include <Eigen/Dense> // for Eigen::
int main(){
// 3*3 matrix A, Bと3*1 vector cの宣言と初期化
Eigen::Matrix3d A;
Eigen::Matrix3d B;
Eigen::Vector3d c;
Eigen::Vector3d d;
Eigen::Matrix3d id = Eigen::Matrix3d::Identity(3, 3); // 単位行列で初期化
A << 2,3,1, 4,-5,-1, 6,-7,1;
B << 1,3,4, -5,-6,3, 6,-1,0;
c << 7, 3, 5;
d << -3, 1, 8;
// matrix A, vecotr c, matrix idの出力
std::cout << A << std::endl;
std::cout << c << std::endl;
std::cout << id << std::endl;
// 要素へのアクセス(0行0列から始まることに注意)
std::cout << A(2, 2) << std::endl; // matrix Aの(2,2)成分
std::cout << A.row(0) << std::endl; // matrix Aの0行目
std::cout << B.col(1) << std::endl; // matrix Bの1列目
std::cout << c.segment(0,2) << std::endl; //vector cの0要素目から1要素まで
// 演算
std::cout << A + B << std::endl; // matrix Aとmatrix Bの和
std::cout << A - B << std::endl; // matrix Aとmatrix Bの差
std::cout << 2 * A << std::endl; // matrix Aの各要素を2倍
std::cout << A * B << std::endl; // matrix Aとmatrix Bの積
// 転置, 随伴, 逆行列
std::cout << A.transpose() << std::endl; // matrix Aの転置行列
std::cout << B.adjoint() << std::endl; // matrix Bの随伴行列
std::cout << A.inverse() << std::endl; // matrix Aの逆行列
// vectorの内積と外積
std::cout << c.dot(d) << std::endl; // vector cとvector dの内積
std::cout << c.cross(d) << std::endl; //vector cとvector dの外積
// Ax = cを解く
Eigen::ColPivHouseholderQR <Eigen::Matrix3d> dec(A);
Eigen::Vector3d x = dec.solve(c);
std::cout << x << std::endl;
// LU分解でAx = cを解く
Eigen::FullPivLU <Eigen::Matrix3d> lu(A); //Eigen::PartialPivLUでも良い
Eigen::Vector3d y = lu.solve(c);
std::cout << y << std::endl;
// 固有値分解
Eigen::EigenSolver <Eigen::Matrix3d> eigensolver(A);
std::cout << eigensolver.eigenvalues() << std::endl; // 固有値
std::cout << eigensolver.eigenvectors() << std::endl; // 固有ベクトル
}