From 0d42b514db5c5aabc47cb83458edd1777809d781 Mon Sep 17 00:00:00 2001 From: chvmvd Date: Sun, 15 Jan 2023 11:21:23 +0900 Subject: [PATCH] Add gaussian-elimination article --- .../10gaussian-elimination/index.mdx | 626 ++++++++++++++++++ .../gauss_jordan_elimination.ipynb | 49 ++ .../gaussian_elimination.ipynb | 47 ++ .../gaussian_elimination_error.ipynb | 52 ++ .../gaussian_elimination_revised.ipynb | 53 ++ 5 files changed, 827 insertions(+) create mode 100644 docs/02algorithms/10gaussian-elimination/index.mdx create mode 100644 static/gaussian-elimination/gauss_jordan_elimination.ipynb create mode 100644 static/gaussian-elimination/gaussian_elimination.ipynb create mode 100644 static/gaussian-elimination/gaussian_elimination_error.ipynb create mode 100644 static/gaussian-elimination/gaussian_elimination_revised.ipynb diff --git a/docs/02algorithms/10gaussian-elimination/index.mdx b/docs/02algorithms/10gaussian-elimination/index.mdx new file mode 100644 index 000000000..2c1313830 --- /dev/null +++ b/docs/02algorithms/10gaussian-elimination/index.mdx @@ -0,0 +1,626 @@ +--- +sidebar_position: 10 +--- + +import ViewSource from "@site/src/components/ViewSource"; +import Answer from "@site/src/components/Answer"; + +# 連立一次方程式の解法 + +$$ +\gdef\customtextcircled#1{\textcircled{\small#1}} +\gdef\formulanumber#1{\,\cdot\cdot\cdot\cdot\cdot\cdot\,\customtextcircled{#1}} +$$ + +## 連立一次方程式の解き方 + +次のような連立一次方程式は、以下のようにして解くと思います。 + +$$ +\gdef\customtextcircled#1{\textcircled{\small#1}} +\gdef\formulanumber#1{\,\cdot\cdot\cdot\cdot\cdot\cdot\,\customtextcircled{#1}} +\begin{dcases} + \begin{alignat*}{3.5} + & x_1 - {} & 2 & x_2 + {} & 3 & x_3 = 3 & \formulanumber{1} \\ + -& x_1 + {} & 3 & x_2 - {} & 2 & x_3 = 1 & \formulanumber{2} \\ + & x_1 - {} & & x_2 + {} & 6 & x_3 = 11 & \formulanumber{3} + \end{alignat*} +\end{dcases} +$$ + +$ +\gdef\customtextcircled#1{\textcircled{\small#1}} +\gdef\formulanumber#1{\,\cdot\cdot\cdot\cdot\cdot\cdot\,\customtextcircled{#1}} +\customtextcircled{1} + \customtextcircled{2}$ より、 + +$$ +\gdef\customtextcircled#1{\textcircled{\small#1}} +\gdef\formulanumber#1{\,\cdot\cdot\cdot\cdot\cdot\cdot\,\customtextcircled{#1}} +x_2 + x_3 = 4 \formulanumber{4} +$$ + +$ +\gdef\customtextcircled#1{\textcircled{\small#1}} +\gdef\formulanumber#1{\,\cdot\cdot\cdot\cdot\cdot\cdot\,\customtextcircled{#1}} +\customtextcircled{2} + \customtextcircled{3}$ より、 + +$$ +\gdef\customtextcircled#1{\textcircled{\small#1}} +\gdef\formulanumber#1{\,\cdot\cdot\cdot\cdot\cdot\cdot\,\customtextcircled{#1}} +\begin{align*} + 2x_2 + 4x_3 &= 12 \\ + \therefore x_2 + 2x_3 &= 6 \formulanumber{5} +\end{align*} +$$ + +$ +\gdef\customtextcircled#1{\textcircled{\small#1}} +\gdef\formulanumber#1{\,\cdot\cdot\cdot\cdot\cdot\cdot\,\customtextcircled{#1}} +\customtextcircled{4},\customtextcircled{5}$ より、 + +$$ +\begin{dcases} + x_2 = 2 \\ + x_3 = 2 +\end{dcases} +$$ + +よって、 + +$$ +\begin{dcases} + x_1 = 1 \\ + x_2 = 2 \\ + x_3 = 2 +\end{dcases} +$$ + +しかし、このようなアルゴリズムでプログラムを作るのは、難しそうです。 + +## Gauss の消去法 + +Gauss の消去法を使えば、システマティックに連立一次方程式を解くことができます。掃き出し法とも呼ばれます。 + +Gauss の消去法は、前進消去と後退代入の二段階から成ります。 + +簡単に Gauss の消去法を説明しておきます。 + +次のように $n$ 個の未知数 $x_1, x_2, x_3, \dots , x_n$ に対して、$m$ 個の方程式を考えます。 + +$$ +\begin{dcases} + \begin{alignat*}{5} + a_{1, 1} & x_1 + {} & a_{1, 2} & x_2 + {} & a_{1, 3} & x_3 + \dots + {} & a_{1, n} & x_n = c_1 \\ + a_{2, 1} & x_1 + {} & a_{2, 2} & x_2 + {} & a_{2, 3} & x_3 + \dots + {} & a_{2, n} & x_n = c_2 \\ + a_{3, 1} & x_1 + {} & a_{3, 2} & x_2 + {} & a_{3, 3} & x_3 + \dots + {} & a_{3, n} & x_n = c_3 \\ + & \cdots\cdot\cdot \\ + a_{m, 1} & x_1 + {} & a_{m, 2} & x_2 + {} & a_{m, 3} & x_3 + \dots + {} & a_{m, n} & x_n = c_m + \end{alignat*} +\end{dcases} +$$ + +この方程式系に対して、以下を行ったものは元の方程式系と同値です。 + +- 二つの方程式を入れ替える +- ある方程式に 0 でないスカラーを掛ける +- ある方程式に他の方程式のスカラー倍を掛ける + +ここで、この方程式系の拡大係数行列を考えます。 + +$$ +\tilde{A} = (A|\bm{c}) = +\left( +\begin{array}{ccccc|c} + a_{1, 1} & a_{1, 2} & a_{1, 3} & \dots & a_{1, n} & c_1 \\ + a_{2, 1} & a_{2, 2} & a_{2, 3} & \dots & a_{2, n} & c_2 \\ + a_{3, 1} & a_{3, 2} & a_{3, 3} & \dots & a_{3, n} & c_3 \\ + \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ + a_{m, 1} & a_{m, 2} & a_{m, 3} & \dots & a_{m, n} & c_m +\end{array} +\right) +$$ + +基本行列を次のように定義します。 + +- $i$ 行と $j$ 行を入れ替える行列を $P_{i, j}$ +- $i$ 行を $\lambda$ 倍する行列を $Q_{i, \lambda}$ +- $i$ 行に $j$ 行の $\lambda$ 倍を加える行列を $R_{i, j, \lambda}$ + + + +拡大係数行列に対して、基本行列を左から掛けて基本変形を繰り返すと、行階段行列 $\tilde{B}$ が作れます。 + +$$ +\tilde{B} = (B|\bm{d}) = +\left( +\begin{array}{cccccccc|c} + 0 & 1 & * & * & * & * & \dots & * & * \\ + 0 & 0 & 0 & 1 & * & * & \dots & * & * \\ + 0 & 0 & 0 & 0 & 0 & 1 & \dots & * & * \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ + 0 & 0 & 0 & 0 & 0 & 0 & \dots & 1 & * \\ + 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 & * \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ + 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 & 0 +\end{array} +\right) +$$ + +これから、作られる方程式系は次のようになりこれははじめの方程式系と同値です。 + +$$ +\begin{dcases} + \begin{alignat*}{4} + & x_{j_1} + \dots + {} & b_{1, j_2} & x_{j_2} + \dots + {} & b_{1, j_3} & x_{j_3} + \dots + {} & b_{1, n} & x_n = d_1 \\ + 0 & x_{j_1} + \dots + {} & & x_{j_2} + \dots + {} & b_{2, j_3} & x_{j_3} + \dots + {} & b_{2, n} & x_n = d_2 \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & & x_{j_3} + \dots + {} & b_{3, n} & x_n = d_3 \\ + & \dots \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & & x_n = d_l \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & 0 & x_n = d_{l + 1} \\ + & \dots \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & 0 & x_n = d_m \\ + \end{alignat*} +\end{dcases} +$$ + +これで $x_n$ から順番に求めていくことで連立方程式を解くことができます。 + +実際に連立方程式を解いてみましょう。 +Gauss の消去法を次の方程式系について行うと、以下のようになります。 + +$$ +\begin{dcases} + \begin{alignat*}{3.5} + & x_1 - {} & 2 & x_2 + {} & 3 & x_3 = 3 \\ + -& x_1 + {} & 3 & x_2 - {} & 2 & x_3 = 1 \\ + & x_1 - {} & & x_2 + {} & 6 & x_3 = 11 + \end{alignat*} +\end{dcases} +$$ + +まずは、前進消去を行います。 + +$$ +\begin{alignat*}{2} + \tilde{A} + & = & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + -1 & 3 & -2 & 1 \\ + 1 & -1 & 6 & 11 + \end{array} + \right) \\ + & \overset{R_{2, 1, 1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + 0 & 1 & 1 & 4 \\ + 1 & -1 & 6 & 11 + \end{array} + \right) \\ + & \overset{R_{3, 1, -1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + 0 & 1 & 1 & 4 \\ + 0 & 1 & 3 & 8 + \end{array} + \right) \\ + & \overset{R_{3, 2, -1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + 0 & 1 & 1 & 4 \\ + 0 & 0 & 2 & 4 + \end{array} + \right) \\ + & \overset{Q_{3, \frac{1}{2}}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + 0 & 1 & 1 & 4 \\ + 0 & 0 & 1 & 2 + \end{array} + \right) +\end{alignat*} +$$ + +次に、連立方程式に戻して後退代入を行っていきます。 + +$$ +\begin{dcases} + \begin{alignat*}{3} + x_1 - 2 & x_2 + {} & 3 & x_3 & {} = {} & 3 \\ + & x_2 + {} & & x_3 & {} = {} & 4 \\ + & & & x_3 & {} = {} & 2 + \end{alignat*} +\end{dcases} \\ +\therefore +\begin{dcases} + \begin{alignat*}{3} + x_1 - 2 & x_2 & & & {} = {} & 3 - 3\times 2 = -3 \\ + & x_2 & & & {} = {} & 4 - 2 = 2 \\ + & & & x_3 & {} = {} & 2 + \end{alignat*} +\end{dcases} \\ +\therefore +\begin{dcases} + \begin{alignat*}{3} + x_1 & & & & {} = {} & -3 + 2\times 2 = 1 \\ + & x_2 & & & {} = {} & 2 \\ + & & & x_3 & {} = {} & 2 + \end{alignat*} +\end{dcases} \\ +\therefore +\begin{dcases} + \begin{align*} + x_1 &= 1\\ + x_2 &= 2\\ + x_3 &= 2 + \end{align*} +\end{dcases} +$$ + +Gauss の消去法を使えば、このようにシステマティックに連立方程式を解けます。 +これなら、プログラムを作るのは簡単そうです。 + +:::info + +次のように、前進消去の段階で行簡約行列を作れば、後退代入を行う必要がなくなります。これは、Gauss-Jordan の消去法と呼ばれます。 +しかし、実は先程のように行簡約行列まで計算しないで途中で止める Gauss の消去法の方が計算量が少なくなるので、Gauss の消去法を使って説明をしました。 + +行簡約行列は次のようになります。 + +$$ +\tilde{B} = +\left( +\begin{array}{cccccccc|c} + 0 & 1 & * & 0 & * & 0 & \dots & * & * \\ + 0 & 0 & 0 & 1 & * & 0 & \dots & * & * \\ + 0 & 0 & 0 & 0 & 0 & 1 & \dots & * & * \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ + 0 & 0 & 0 & 0 & 0 & 0 & \dots & 1 & * \\ + 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 & * \\ + \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ + 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 & 0 +\end{array} +\right) +$$ + +これから、作られる方程式系は次のようになりこれは元の方程式系と同値です。 + +$$ +\begin{dcases} + \begin{alignat*}{4} + & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & 0 & x_n = d_1 \\ + 0 & x_{j_1} + \dots + {} & & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & 0 & x_n = d_2 \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & & x_{j_3} + \dots + {} & 0 & x_n = d_3 \\ + & \dots \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & & x_n = d_l \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & 0 & x_n = d_{l + 1} \\ + & \dots \\ + 0 & x_{j_1} + \dots + {} & 0 & x_{j_2} + \dots + {} & 0 & x_{j_3} + \dots + {} & 0 & x_n = d_m \\ + \end{alignat*} +\end{dcases} +$$ + +これで連立方程式を解くことができました。 + +実際に連立方程式を解いてみます。 +これを次の方程式系について行うと、以下のようになります。 + +$$ +\begin{dcases} + \begin{alignat*}{3.5} + & x_1 - {} & 2 & x_2 + {} & 3 & x_3 = 3 \\ + -& x_1 + {} & 3 & x_2 - {} & 2 & x_3 = 1 \\ + & x_1 - {} & & x_2 + {} & 6 & x_3 = 11 + \end{alignat*} +\end{dcases} +$$ + +$$ +\begin{alignat*}{5} + \tilde{A} + & = & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + -1 & 3 & -2 & 1 \\ + 1 & -1 & 6 & 11 + \end{array} + \right) \\ + & \overset{R_{2, 1, 1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + 0 & 1 & 1 & 4 \\ + 1 & -1 & 6 & 11 + \end{array} + \right) \\ + & \overset{R_{3, 1, -1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -2 & 3 & 3 \\ + 0 & 1 & 1 & 4 \\ + 0 & 1 & 3 & 8 + \end{array} + \right) \\ + & \overset{R_{1, 2, 2}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & 0 & 5 & 11 \\ + 0 & 1 & 1 & 4 \\ + 0 & 1 & 3 & 8 + \end{array} + \right) \\ + & \overset{R_{3, 2, -1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & 0 & 5 & 11 \\ + 0 & 1 & 1 & 4 \\ + 0 & 0 & 2 & 4 + \end{array} + \right) \\ + & \overset{Q_{3, \frac{1}{2}}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & 0 & 5 & 11 \\ + 0 & 1 & 1 & 4 \\ + 0 & 0 & 1 & 2 + \end{array} + \right) \\ + & \overset{R_{1, 3, -5}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & 0 & 0 & 1 \\ + 0 & 1 & 1 & 4 \\ + 0 & 0 & 1 & 2 + \end{array} + \right) \\ + & \overset{R_{2, 3, -1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & 0 & 0 & 1 \\ + 0 & 1 & 0 & 2 \\ + 0 & 0 & 1 & 2 + \end{array} + \right) \\ +\end{alignat*} +$$ + +$$ +\begin{dcases} + \begin{alignat*}{4} + & x_1 + {} & 0 & x_2 + {} & 0 & x_3 & & = 1 \\ + 0 & x_1 + {} & & x_2 + {} & 0 & x_3 & & = 2 \\ + 0 & x_1 + {} & 0 & x_2 + {} & & x_3 & & = 2 + \end{alignat*} +\end{dcases} \\ +\therefore +\begin{dcases} + \begin{align*} + x_1 &= 1 \\ + x_2 &= 2 \\ + x_3 &= 2 + \end{align*} +\end{dcases} +$$ + +::: + +## Gauss の消去法を使ったプログラム + +Gauss の消去法を使って連立一次方程式を解くプログラムを作ってみましょう。概要を説明します。 + +まずは、前進消去を行います。 + +$i$ 行、$i$ 列が pivot となるので、$i$ 行目を pivot の値で割って pivot を 1 にします。 +その後、その行で $i+1$ 行目以降を掃き出します。 + +次は、後退代入を行います。 + +$x_n$ は $d_n$ になります。求まった $x_n$ をそれよりも上の式に代入して、$b_{j, n - 1}x_{n - 1}$ を右辺に移動させ、$d_{n - 1}$ の値を更新します。これを繰り返すと、後退代入ができます。 + +プログラムは次のようになります。計算量は、$O(n^3)$ です。 + + + +`reversed` 関数は配列の要素を反転します。 + +## 部分ピボット選択 + +次のような連立方程式を解こうとすると、次のようなエラーが出てしまいます。 + +$$ +\begin{dcases} + \begin{alignat*}{3.5} + & & -2 & x_2 + {} & 3 & x_3 = 2 \\ + -& x_1 + {} & 3 & x_2 - {} & 2 & x_3 = 1 \\ + & x_1 - {} & & x_2 + {} & 6 & x_3 = 11 + \end{alignat*} +\end{dcases} +$$ + + + +これは、前進消去の際に 0 で割る操作ができてしまったためです。 +これを解決するために、部分ピボット選択を行います。誤差を小さくする役割もあります。 + +部分ピボット選択は、pivot 列の pivot 行以降の行で絶対値が最大になる行を pivot に使うように変形することです。 + +先程の連立方程式なら、次のようになります。 + +$$ +\begin{alignat*}{2} + & & + & \left( + \begin{array}{ccc|c} + 0 & -2 & 3 & 2 \\ + -1 & 3 & -2 & 1 \\ + 1 & -1 & 6 & 11 + \end{array} + \right) \\ + & \overset{P_{1, 2}}{\to} & + & \left( + \begin{array}{ccc|c} + -1 & 3 & -2 & 1 \\ + 0 & -2 & 3 & 2 \\ + 1 & -1 & 6 & 11 + \end{array} + \right) \\ + & \overset{Q_{2, -1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -3 & 2 & -1 \\ + 0 & -2 & 3 & 2 \\ + 1 & -1 & 6 & 11 + \end{array} + \right) \\ + & \overset{R_{3, 1, -1}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -3 & 2 & -1 \\ + 0 & -2 & 3 & 2 \\ + 0 & 2 & 4 & 12 + \end{array} + \right) \\ + & \overset{Q_{2, -\frac{1}{2}}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -3 & 2 & -1 \\ + 0 & 1 & -\frac{3}{2} & -1 \\ + 0 & 2 & 4 & 12 + \end{array} + \right) \\ + & \overset{R_{3, 2, -2}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -3 & 2 & -1 \\ + 0 & 1 & -\frac{3}{2} & -1 \\ + 0 & 0 & 7 & 14 + \end{array} + \right) \\ + & \overset{Q_{3, \frac{1}{7}}}{\to} & + & \left( + \begin{array}{ccc|c} + 1 & -3 & 2 & -1 \\ + 0 & 1 & -\frac{3}{2} & -1 \\ + 0 & 0 & 1 & 2 + \end{array} + \right) \\ +\end{alignat*} +$$ + +部分ピボット選択を使うと、次のようになります。 + + + +## 練習問題 + +Gauss-Jordan の消去法を使って行簡約行列を作ることで、連立方程式を解くプログラムを作ってください。 + + + + + + diff --git a/static/gaussian-elimination/gauss_jordan_elimination.ipynb b/static/gaussian-elimination/gauss_jordan_elimination.ipynb new file mode 100644 index 000000000..e60d6cb29 --- /dev/null +++ b/static/gaussian-elimination/gauss_jordan_elimination.ipynb @@ -0,0 +1,49 @@ +{ + "nbformat": 4, + "nbformat_minor": 2, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "def gausu_jordan_elimination(a):\n", + " # 前進消去\n", + " for i in range(len(a)):\n", + " # 部分ピボット選択\n", + " for j in range(i + 1, len(a)):\n", + " if abs(a[i][i]) < abs(a[j][i]):\n", + " tmp = a[i]\n", + " a[i] = a[j]\n", + " a[j] = tmp\n", + " # pivot倍で行を割る\n", + " pivot = a[i][i]\n", + " for j in range(i, len(a[i])):\n", + " a[i][j] /= pivot\n", + " # 掃き出す\n", + " for j in range(len(a)):\n", + " if j != i:\n", + " factor = a[j][i]\n", + " for k in range(i, len(a[i])):\n", + " a[j][k] -= factor * a[i][k]\n", + " x = []\n", + " for i in range(len(a)):\n", + " x.append(a[i][len(a[i]) - 1])\n", + " return x\n", + "\n", + "\n", + "print(gausu_jordan_elimination([[0, -2, 3, 2], [-1, 3, -2, 1], [1, -1, 6, 11]]))" + ], + "cell_type": "code", + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[1.0, 2.0, 2.0]\n" + ] + } + ], + "execution_count": null + } + ] +} diff --git a/static/gaussian-elimination/gaussian_elimination.ipynb b/static/gaussian-elimination/gaussian_elimination.ipynb new file mode 100644 index 000000000..2a009554b --- /dev/null +++ b/static/gaussian-elimination/gaussian_elimination.ipynb @@ -0,0 +1,47 @@ +{ + "nbformat": 4, + "nbformat_minor": 2, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "def gaussian_elimination(a):\n", + " # 前進消去\n", + " for i in range(len(a)):\n", + " # pivot倍で行を割る\n", + " pivot = a[i][i]\n", + " for j in range(i, len(a[i])):\n", + " a[i][j] /= pivot\n", + "\n", + " # i+1行目以降を掃き出す\n", + " for j in range(i + 1, len(a)):\n", + " factor = a[j][i]\n", + " for k in range(i, len(a[i])):\n", + " a[j][k] -= factor * a[i][k]\n", + " # 後退代入\n", + " x = []\n", + " for i in reversed(range(len(a))):\n", + " tmp = a[i][len(a[i]) - 1] / a[i][i]\n", + " for j in range(i):\n", + " a[j][len(a[i]) - 1] -= a[j][i] * tmp\n", + " x = [tmp] + x\n", + " return x\n", + "\n", + "\n", + "print(gaussian_elimination([[1, -2, 3, 3], [-1, 3, -2, 1], [1, -1, 6, 11]]))" + ], + "cell_type": "code", + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[1.0, 2.0, 2.0]\n" + ] + } + ], + "execution_count": null + } + ] +} diff --git a/static/gaussian-elimination/gaussian_elimination_error.ipynb b/static/gaussian-elimination/gaussian_elimination_error.ipynb new file mode 100644 index 000000000..4d7d0a35b --- /dev/null +++ b/static/gaussian-elimination/gaussian_elimination_error.ipynb @@ -0,0 +1,52 @@ +{ + "nbformat": 4, + "nbformat_minor": 2, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "def gaussian_elimination(a):\n", + " # 前進消去\n", + " for i in range(len(a)):\n", + " # pivot倍で行を割る\n", + " pivot = a[i][i]\n", + " for j in range(i, len(a[i])):\n", + " a[i][j] /= pivot\n", + "\n", + " # i+1行目以降を掃き出す\n", + " for j in range(i + 1, len(a)):\n", + " factor = a[j][i]\n", + " for k in range(i, len(a[i])):\n", + " a[j][k] -= factor * a[i][k]\n", + " # 後退代入\n", + " x = []\n", + " for i in reversed(range(len(a))):\n", + " tmp = a[i][len(a[i]) - 1] / a[i][i]\n", + " for j in range(i):\n", + " a[j][len(a[i]) - 1] -= a[j][i] * tmp\n", + " x = [tmp] + x\n", + " return x\n", + "\n", + "\n", + "print(gaussian_elimination([[0, -2, 3, 2], [-1, 3, -2, 1], [1, -1, 6, 11]]))" + ], + "cell_type": "code", + "outputs": [ + { + "output_type": "error", + "ename": "ZeroDivisionError", + "evalue": "division by zero", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mZeroDivisionError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn [1], line 23\u001b[0m\n\u001b[1;32m 20\u001b[0m x\u001b[39m=\u001b[39m[tmp]\u001b[39m+\u001b[39mx\n\u001b[1;32m 21\u001b[0m \u001b[39mreturn\u001b[39;00m x\n\u001b[0;32m---> 23\u001b[0m \u001b[39mprint\u001b[39m(gaussian_elimination([[\u001b[39m0\u001b[39m,\u001b[39m-\u001b[39m\u001b[39m2\u001b[39m,\u001b[39m3\u001b[39m,\u001b[39m2\u001b[39m],[\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m,\u001b[39m3\u001b[39m,\u001b[39m-\u001b[39m\u001b[39m2\u001b[39m,\u001b[39m1\u001b[39m],[\u001b[39m1\u001b[39m,\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m,\u001b[39m6\u001b[39m,\u001b[39m11\u001b[39m]]))\n", + "Cell \u001b[0;32mIn [1], line 7\u001b[0m, in \u001b[0;36mgaussian_elimination\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m 5\u001b[0m pivot\u001b[39m=\u001b[39ma[i][i]\n\u001b[1;32m 6\u001b[0m \u001b[39mfor\u001b[39;00m j \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(i,\u001b[39mlen\u001b[39m(a[i])):\n\u001b[0;32m----> 7\u001b[0m a[i][j]\u001b[39m/\u001b[39m\u001b[39m=\u001b[39mpivot\n\u001b[1;32m 9\u001b[0m \u001b[39m# i+1行目以降を掃き出す\u001b[39;00m\n\u001b[1;32m 10\u001b[0m \u001b[39mfor\u001b[39;00m j \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(i\u001b[39m+\u001b[39m\u001b[39m1\u001b[39m,\u001b[39mlen\u001b[39m(a)):\n", + "\u001b[0;31mZeroDivisionError\u001b[0m: division by zero" + ] + } + ], + "execution_count": null + } + ] +} diff --git a/static/gaussian-elimination/gaussian_elimination_revised.ipynb b/static/gaussian-elimination/gaussian_elimination_revised.ipynb new file mode 100644 index 000000000..4e394b839 --- /dev/null +++ b/static/gaussian-elimination/gaussian_elimination_revised.ipynb @@ -0,0 +1,53 @@ +{ + "nbformat": 4, + "nbformat_minor": 2, + "metadata": {}, + "cells": [ + { + "metadata": {}, + "source": [ + "def gaussian_elimination_revised(a):\n", + " # 前進消去\n", + " for i in range(len(a)):\n", + " # 部分ピボット選択\n", + " for j in range(i + 1, len(a)):\n", + " if abs(a[i][i]) < abs(a[j][i]):\n", + " tmp = a[i]\n", + " a[i] = a[j]\n", + " a[j] = tmp\n", + " # pivot倍で行を割る\n", + " pivot = a[i][i]\n", + " for j in range(i, len(a[i])):\n", + " a[i][j] /= pivot\n", + "\n", + " # i+1行目以降を掃き出す\n", + " for j in range(i + 1, len(a)):\n", + " factor = a[j][i]\n", + " for k in range(i, len(a[i])):\n", + " a[j][k] -= factor * a[i][k]\n", + " # 後退代入\n", + " x = []\n", + " for i in reversed(range(len(a))):\n", + " tmp = a[i][len(a[i]) - 1] / a[i][i]\n", + " for j in range(i):\n", + " a[j][len(a[i]) - 1] -= a[j][i] * tmp\n", + " x = [tmp] + x\n", + " return x\n", + "\n", + "\n", + "print(gaussian_elimination_revised([[0, -2, 3, 2], [-1, 3, -2, 1], [1, -1, 6, 11]]))" + ], + "cell_type": "code", + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[1.0, 2.0, 2.0]\n" + ] + } + ], + "execution_count": null + } + ] +}