はじめに

プログラミングで二元一次方程式を解く方法を紹介します．
泥臭く数式変形する解法とオシャレに行列を使用する解法を紹介します．
私自身プログラミング歴が浅いのもあり，プログラムによくない書き方をしていることも考えられます．何かお気付きでしたらご指摘いただけると嬉しいです．

解く数式

ここで解法を求める数式は以下のものとする．

\left\{
\begin{array}{ll}
ax + by = c \quad (1)\\
dx + ey = f \quad (2)
\end{array}
\right.

数式変形による解法

ここでは，泥臭く式変形をして解の式にする．式変形に手順は以下．
まずは，$(1)\times e - (2)\times b$より

(ae-bd)x=(ce-bf)\\
x = \frac{ce-bf}{ae-bd}\quad (3)

$(1)$より，

y=\frac{c}{b}-\frac{a}{b}x\quad (4)

$(3)$および$(4)$を用いて方程式を解くプログラムが以下である．ここでは，変数の例として$a=1, b=2, c=1, d=3, e=4, f=2$を用いた．この値における解は$x=0, y=0.5$である．

formula_deformation_solve

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

vector<double> solve(const vector<vector<double>>& m, const vector<double>& v);
void view(const vector<double>& v);

int main(){
  vector<vector<double>> A(2, vector<double>(2, 0.0));
  A[0][0] = 1;
  A[0][1] = 2;
  A[1][0] = 3;
  A[1][1] = 4;

  vector<double> B(2, 0.0);
  B[0] = 1;
  B[1] = 2;

  auto X = solve(A, B);
  view(X);
}

vector<double> solve(const vector<vector<double>>& m, const vector<double>& v){
  vector<double> sol(2);
  sol[0] = (v[0]*m[1][1]-m[0][1]*v[1])/(m[0][0]*m[1][1]-m[0][1]*m[1][0]);
  sol[1] = (v[0]-m[0][0]*sol[0])/m[0][1];
  return sol;
}

void view(const vector<double>& v){
  for(const auto& e : v)
    cout << e << " ";
  cout << endl;
}

実行結果が以下である．成功した．

output

-0 0.5

行列計算による解法

数式変形による方法でも上記のように解を求めることはできるのだが，それでは少々カッコが悪い．そこで，行列計算による解法を考えた．考え方は以下である．$(1)$および$(2)$における各変数によって構成される行列を以下のように定義する．

\boldsymbol{A}=
\begin{bmatrix}
a & b \\
d & e 
\end{bmatrix}
\boldsymbol{B}=
\begin{bmatrix}
c \\
f 
\end{bmatrix}
\boldsymbol{X}=
\begin{bmatrix}
x \\
y 
\end{bmatrix}

こうすることで，$(1)$および$(2)$式が以下のように表される．

\boldsymbol{A}\cdot \boldsymbol{X}=\boldsymbol{B}

ここで，両辺に左から$\boldsymbol{A}^{-1}$を乗じると，

\boldsymbol{A}^{-1}\cdot\boldsymbol{A}\cdot \boldsymbol{X}=
\boldsymbol{A}^{-1}\cdot\boldsymbol{B}
\\
\boldsymbol{X}=\boldsymbol{A}^{-1}\cdot\boldsymbol{B}\quad (5)

$(5)$を用いて方程式を解くプログラムが以下である．ここで，逆行列に関しては掃き出し法など良いアルゴリズムがあるのだが，今回は$2\times 2$逆行列しか求めないため，直接要素を算出する方法を用いた．先ほどと同様，変数の例として$a=1, b=2, c=1, d=3, e=4, f=2$を用いた．この値における解は$x=0, y=0.5$である．

matrix_solve

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

vector<vector<double>> inverse(const vector<vector<double>>& m);
vector<double> dot(const vector<vector<double>>& m, const vector<double>& v);
void view(const vector<double>& v);

int main(){
  vector<vector<double>> A(2, vector<double>(2, 0.0));
  A[0][0] = 1;
  A[0][1] = 2;
  A[1][0] = 3;
  A[1][1] = 4;

  vector<double> B(2, 0.0);
  B[0] = 1;
  B[1] = 2;

  auto A_inv = inverse(A);
  auto X = dot(A_inv, B);
  view(X);
}

vector<vector<double>> inverse(const vector<vector<double>>& m){
  vector<vector<double>> m_inv(m.size(), vector<double>(m.front().size()));
  const double det = 1.0/(m[0][0]*m[1][1]-m[0][1]*m[1][0]);
  m_inv[0][0] = det * m[1][1];
  m_inv[0][1] = - det * m[0][1];
  m_inv[1][0] = - det * m[1][0];
  m_inv[1][1] = det * m[0][0];
  return m_inv;
}

vector<double> dot(const vector<vector<double>>& m, const vector<double>& v){
  vector<double> product(m.size(), 0.0);
  if(m.size()!=v.size()){
    cout << "matrix size and vector size are different." << endl;
    return product;
  }
  for(size_t i=0; i<m.size(); i++)
    for(size_t j=0; j<m[i].size(); j++)
      product[i] += m[i][j] * v[j];
  return product;
}

void view(const vector<double>& v){
  for(const auto& e : v)
    cout << e << " ";
  cout << endl;
}

実行結果が以下である．成功した．

output

-0 0.5

任意の変数で解を求める

上記の例では，変数$a, b, c, d, e, f$を変更する度にコンパイルする必要があるため，不便である．
実行時に標準入力から変数を取得してそれを行列によって解くプログラムが以下である．

arbitrary_variable

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

vector<vector<double>> inverse(const vector<vector<double>>& m);
vector<double> dot(const vector<vector<double>>& m, const vector<double>& v);
void view(const vector<double>& v);

int main(){
  vector<vector<double>> A(2, vector<double>(2, 0.0));
  vector<double> B(2, 0.0);
  cout << "Plese input first equation variable." << endl;
  cin >> A[0][0] >> A[0][1] >> B[0];
  cout << "Plese input second equation variable." << endl;
  cin >> A[1][0] >> A[1][1] >> B[1];

  auto A_inv = inverse(A);
  auto X = dot(A_inv, B);
  cout << "solutions is : " << endl;
  view(X);
}

vector<vector<double>> inverse(const vector<vector<double>>& m){
  vector<vector<double>> m_inv(m.size(), vector<double>(m.front().size()));
  const double det = 1.0/(m[0][0]*m[1][1]-m[0][1]*m[1][0]);
  m_inv[0][0] = det * m[1][1];
  m_inv[0][1] = - det * m[0][1];
  m_inv[1][0] = - det * m[1][0];
  m_inv[1][1] = det * m[0][0];
  return m_inv;
}

vector<double> dot(const vector<vector<double>>& m, const vector<double>& v){
  vector<double> product(m.size(), 0.0);
  if(m.size()!=v.size()){
    cout << "matrix size and vector size are different." << endl;
    return product;
  }
  for(size_t i=0; i<m.size(); i++)
    for(size_t j=0; j<m[i].size(); j++)
      product[i] += m[i][j] * v[j];
  return product;
}

void view(const vector<double>& v){
  for(const auto& e : v)
    cout << e << " ";
  cout << endl;
}

本プログラムでの各値入力順は以下である．値の区切りはスペースまたはエンターとする．

input_order

a b c
d e f

実際の入出力例が以下で，方程式を正しく解くことができた．

sample_input

Plese input first equation variable.
1 2 3
Plese input second equation variable.
4 5 6

sample_output

solutions is : 
-1 2

コメントで教えていただいた方法（行列)

コメントの書き方でmatrixを使った方法が以下である．

matrix_solve_from_comment

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

class matrix {
 public:
  matrix() = default;
  matrix(std::size_t rows, std::size_t cols) : data_(rows*cols), rows_{rows} {}

  size_t size() const noexcept { return data_.size(); }
  size_t rows() const noexcept { return rows_; }
  size_t cols() const noexcept { return this->size()/rows_; }

  void resize(std::size_t rows, std::size_t cols) {
    data_.resize(rows*cols);
    rows_ = rows;
  }

  auto& operator()(std::size_t i, std::size_t j) {
    return data_.at(rows_*j + i);
  }
  const auto& operator()(std::size_t i, std::size_t j) const {
    return data_.at(rows_*j + i);
  }

 private:
  std::vector<double> data_;
  std::size_t rows_;
};

matrix inverse(const matrix& m);
vector<double> dot(const matrix& m, const vector<double>& v);
void view(const matrix& m);
void view(const vector<double>& v);

int main(){
  matrix A{2, 2};
  cin >> A(0,0) >> A(0,1) >> A(1,0) >> A(1,1);
  vector<double> B(2, 0.0);
  cin >> B[0] >> B[1];

  auto A_inv = inverse(A);
  auto X = dot(A_inv, B);
  view(X);
}

matrix inverse(const matrix& m){
  matrix m_inv{m.rows(), m.cols()};
  const double det = 1.0/(m(0,0)*m(1,1)-m(0,1)*m(1,0));
  m_inv(0,0) = det * m(1,1);
  m_inv(0,1) = - det * m(0,1);
  m_inv(1,0) = - det * m(1,0);
  m_inv(1,1) = det * m(0,0);
  return m_inv;
}

vector<double> dot(const matrix& m, const vector<double>& v){
  vector<double> product(m.cols(), 0.0);
  if(m.cols()!=v.size()){
    cout << "matrix size and vector size are different." << endl;
    return product;
  }

  for(std::size_t i=0; i<m.cols(); i++){
    for(std::size_t j=0; j<m.rows(); j++){
      product[i] += m(i,j) * v[j];
    }
  }
  return product;
}

void view(const matrix& m){
  for(std::size_t i=0; i<m.cols(); i++){
    for(std::size_t j=0; j<m.rows(); j++){
      cout << m(i,j) << " ";
    }
    cout << endl;
  }
  cout << endl;
}


void view(const vector<double>& v){
  for(const auto& e : v)
    cout << e << " ";
  cout << endl;
}

おわりに

今回のように簡単な行列では数式変形を用いた解法でも良いのですが，複雑な式について数式を変形するのは骨が折れるので，行列が大事だということを実感できました．