C++で数値計算 (偏微分方程式)

数値計算の基本をC++で書いていきます. 
[YonohaのGithub](https://github.com/Yonoha/numerical_analysis_c_plus)にも全体のコードとグラフを公開しています. 
 
`boost`や`Eigen`, `gnuplot`を使用しているコードのコンパイルには, `boost`や`Eigen`, `gnuplot`のインストールがそれぞれ必要です. 

変更した方が良い箇所があれば, この記事かTwitterの[@Yonono_01](https://twitter.com/Yonono_01)まで連絡ください.

## 目次

- [偏微分方程式 (Partial Differential Equation)](#偏微分方程式-partial-differential-equation)
	- [差分法 (Difference Method)](#差分法-difference-method)
	- [有限要素法 (Finite Element Method)](#有限要素法-finite-element-method)
		- [1次元Poisson方程式](#1次元poisson方程式)
		- [1次元定常移流拡散方程式](#1次元定常移流拡散方程式)
		- [1次元非定常拡散方程式](#1次元非定常拡散方程式)
		- [1次元非定常移流拡散方程式](#1次元非定常移流拡散方程式)
		- [1次元Burgers方程式](#1次元burgers方程式)
	- [その他の解法](#その他の解法)
- [参考文献](#参考文献)

# 偏微分方程式 (Partial Differential Equation)

線形偏微分方程式は楕円型/, 双曲型, 放物型のいずれかに分類される. 
楕円型偏微分方程式にはPoisson方程式(において$f=0$としたLaplace方程式)やHelmholtz方程式などがある. 
双曲型偏微分方程式には波動方程式や移流方程式などがある.
放物型偏微分方程式には拡散方程式やSchrödinger方程式などがある. 
非線形偏微分方程式にはNavier Stokes方程式やKdV方程式などがある. 

この記事では有限要素法を中心に説明する. 

## 差分法 (Difference Method)

次の$1$次元拡散方程式
```math
\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} ~ (0 < x < 1, 0 < t < \infty)
```
を初期条件
```math
u(x,0) = \sin(\pi x) ~ (0 \leq x \leq 1)
```
と境界条件
```math
u(0,t) = u(1,t) = 0 ~ (0 < t < \infty)
```
の元で解く.
この偏微分方程式をメッシュ数を$I$として, 差分スキームの形に書きなおすと

差分方程式
```math
u_i^{n+1} = u_i^n + \frac{d t}{dx^2} (u_{j+1}^n - 2u_j^n + u_{j-1}^n) ~ (n = 0, 1, \cdots, i = 1, \cdots, I-1)
```
初期条件
```math
u_i^0 = \sin(\pi x_i)
```
境界条件
```math
u_0^n = u_I^n = 0
```
のようになる. 

```cpp
#include <cmath>    // for std::sin
#include <fstream>  // for std::ofstream 
#include <iostream> // for std::cerr, std::endl
#include <vector>   // for std::vector
#include <boost/assert.hpp>                   // for BOOST_ASSERT
#include <boost/math/constants/constants.hpp> // for boost::math::constants::pi

namespace{
    using myvec = std::vector<double>;

    static auto constexpr TEND = 10.0;
    static auto constexpr TLOOP = 10000;
    static auto constexpr DT = TEND / TLOOP;
    static auto constexpr TREP = TLOOP / 1000;

    static auto constexpr MESH = 100;
    static auto constexpr XRANGE = 1.0;
    static auto constexpr DX = XRANGE / MESH;

    // boundary condition
    static auto constexpr BOUND_LEFT = 0.0;
    static auto constexpr BOUND_RIGHT = 0.0;
    
    // coefficient
    static auto constexpr DIFF = 0.01;
    static auto constexpr KAP = DIFF * DT / (DX * DX);
    
    // result output file
    bool output_result(std::vector<myvec> const &uxt);

    // discretization of x
    myvec discrete();

    // initial condition u(x, 0)
    myvec initial(myvec const &x);

    // center differential method
    std::vector<myvec> differencial_method(myvec const &u0);
}

int main(){
    // KAP > 0.5のとき強制終了 
    // Neumannの安定性解析より, KAP <= 0.5でないと正しい結果が得られない
    BOOST_ASSERT(KAP <= 0.5);

    auto x = discrete();

    auto u0 = initial(x);

    auto uxt = differencial_method(u0);

    output_result(uxt);

    if (!output_result(uxt)){
        std::cerr << "error! output file not open" << std::endl;
        return -1;
    }

    return 0;
}

namespace{
    bool output_result(std::vector<myvec> const &uxt){
        std::vector<myvec> u(uxt.size(), (myvec(TREP + 1)));
        
        // output only TREP + 1 (TLOOP = 0, 1000, ... , 10000)
        for (auto i = 0; i <= TREP; i++){
            for (auto j = 0; j <= MESH; j++){
                u[j][i] = uxt[j][1000 * i];
            }
        }

        std::ofstream ofs;
        ofs.open("data_center_diffe_diff.txt");

        if (!ofs){
            return false;
        }

        for (auto j = 0; j <= MESH; j++){
            ofs << static_cast<double>(j) * DX << " ";
            
            for (auto i = 0; i <= TREP; i++){
                ofs << u[j][i] << " ";
            }
            ofs << std::endl;
        }

        return true;
    }

    myvec discrete(){
        myvec x(MESH + 1, 0.0);

        // discrization of xrange
        for (int i = 0; i <= MESH; i++){
            x[i] = i * DX;
        }

        return x;
    }

    myvec initial(myvec const &x){
        myvec u0(MESH + 1, 0.0);

        // initial condition
        for (int i = 0; i<= MESH; i++){
            u0[i] = std::sin(boost::math::constants::pi<double>() * x[i]);
        }

        return u0;
    }

    std::vector<myvec> differencial_method(myvec const &u0){
        std::vector<myvec> uxt(MESH + 1, (myvec(TLOOP + 1)));

        // initial condition
        for (auto j = 0; j <= MESH; j++){
            uxt[j][0] = u0[j];
        }

        for (auto i = 0; i < TLOOP; i++){
            // boundary condition
            uxt[0][i + 1] = BOUND_LEFT;
            uxt[MESH][i + 1] = BOUND_RIGHT;

            // center differential method
            for (auto j = 1; j < MESH; j++){
                uxt[j][i + 1] = uxt[j][i] + KAP * (uxt[j + 1][i] - 2.0 * uxt[j][i] + uxt[j - 1][i]);
            }
        }

        return uxt;
    }
}
```

出力ファイル`data_center_diffe_diff.txt`からデータを読み取り, gnuplotを用いてグラフを描き, `center_diffe_diff.png`に保存する.

```gnuplot
set title '1-dim diffusion equation by center differential method'

set xl 'x'
set yl 'u'
set grid

plot for [i = 2 : 12 : 2] 'data_center_diffe_diff.txt' using 1:i w l notitle

set terminal pngcairo
set output 'center_diffe_diff.png'
replot
```

`center_diffe_diff.png`は次のようになる. 

![center_diffe_diff](https://user-images.githubusercontent.com/75283304/175845164-ab9ed5ce-c1fc-425f-a8e8-56e18ef1bdc8.png)

## 有限要素法 (Finite Element Method)

有限要素法は微分方程式の弱形式をもとにした近似解法である.

Dirichlet条件, 左Neumann条件, 右Neumann条件のいずれにも対応するようにコードを書く. 

### 1次元Poisson方程式

$f$を既知関数とし, $u$を求める未知関数とする. 
このとき, 次の偏微分方程式
```math
\Delta u + f = 0
```
をPoisson方程式という. 
$1$次元においては
```math
\frac{d^2 u}{dx^2} + f = 0 ~ (x \in [0, L])
```
となる. 
今回は
```math
f = 1
```
の場合を考える. 

```cpp
#include <fstream>  // for std::ofstream
#include <iostream> // for std::endl
#include <vector>   // fos std::vector
#include <boost/assert.hpp> // for BOOST_ASSERT

namespace fem{
	using myvec = std::vector<double>;

    class FEM final{
    	private:
        	enum class boundary_condi_type{
            	DIRICLET = 0,
            	LEFT_NEUMANN = 1,
            	RIGHT_NEUMANN = 2
        	};
        
        	static auto constexpr BCT = boundary_condi_type::LEFT_NEUMANN;

			// boundary conditions
        	static auto constexpr D0 = 0.0; // left Dirichlet
    		static auto constexpr D1 = 1.0; // right Direchlet
			static auto constexpr N0 = -1.0; // left Neumann
	    	static auto constexpr N1 = 0.5; // right Neumann

	    	static auto constexpr ELEMENT = 100;
        	static auto constexpr LENGTH = 1.0;
	    	static auto constexpr NODE = ELEMENT + 1;
        	static auto constexpr DX = LENGTH / ELEMENT;
	
    	private:
        	myvec bound_;
        	myvec diag_;
			myvec f_;
	    	myvec left_;
	    	myvec right_;

    	public:
        	FEM()
            	: bound_(NODE, 0.0), diag_(NODE, 0.0),f_(NODE), left_(NODE, 0.0), right_(NODE, 0.0) 
			{
				for (auto && elem : f_){
					elem = 1.0; // f = 1
				}

			}

        	~FEM() = default;

			// コピーコントラクタでcopy禁止
        	FEM(FEM const &dummy) = delete;

			// operator=()でもcopy禁止
        	FEM & operator=(FEM const &dummy) = delete;

			// boundary condition
        	void boundary();

			// make stiffness matrix
	    	void mat();;

			// TDMA method
	    	myvec tdma() const;

			// result output file
			bool output_file(myvec const &x);
    };
}

int main(){
    fem::FEM fem_obj;

	fem_obj.mat();
	fem_obj.boundary();
	auto const x = fem_obj.tdma();

	if (!fem_obj.output_file(x)){
        std::cerr << "output file not open" << std::endl;
        return -1;
    }

	return 0;
}

namespace fem{
	bool FEM::output_file(myvec const &x){
		std::ofstream ofs("data_Poisson.txt");
    	if (!ofs) {
        	return false;
    	}

		auto const size = static_cast<int>(x.size());

		for (auto i = 0; i < size; i++){
			ofs << fem::FEM::DX * static_cast<double>(i) << " " << x[i] << std::endl;
		}

		return true;
	}

	myvec FEM::tdma() const{
		myvec p(NODE, 0.0);
		myvec q(NODE, 0.0);

		p[0] = -right_[0] / diag_[0];
		q[0] = bound_[0] / diag_[0];

		for (auto i = 1; i < NODE; i++){
		    p[i] = -right_[i] / (diag_[i] + left_[i] * p[i - 1]);
		    q[i] = (bound_[i] - left_[i] * q[i - 1]) / (diag_[i] + left_[i] * p[i - 1]);
	    }

		myvec x(NODE, 0.0);
		x[NODE - 1] = q[NODE - 1];

		for (auto j = NODE - 2; j >= 0; j--){
			x[j] = p[j] * x[j + 1] + q[j];
		}
		
        return x;
 	}

    void FEM::boundary(){
		switch (BCT){
            case boundary_condi_type::DIRICLET:
            	diag_[0] = 1.0;
				diag_[NODE - 1] = 1.0;
			    left_[NODE - 1] = 0.0;
        	    right_[0] = 0.0;
				// boundary conditions
				bound_[0] = D0;
        	    bound_[NODE - 1] = D1;
                break;
		
            case boundary_condi_type::LEFT_NEUMANN:
        	    diag_[NODE - 1] = 1.0;
				left_[NODE - 1] = 0.0;
				// boundary conditions
			    bound_[0] -= N0; // left Neumann
				bound_[NODE - 1] = D1;
                break;
			
            case boundary_condi_type::RIGHT_NEUMANN:
			    diag_[0] = 1.0;
        	    right_[0] = 0.0;
				// boundary conditions
        	    bound_[0] = D0;
			    bound_[NODE - 1] += N1; // right Neumann
                break;
            
            default:
                BOOST_ASSERT(!"error! wrong boundary condition");
                break;
        }
	}

	void FEM::mat(){
	    for (auto i = 0; i < ELEMENT; i++){
			//boundary
		    bound_[i] += (2.0 * f_[i] + 1.0 * f_[i + 1]) * DX / 6.0;
		    bound_[i + 1] += (1.0 * f_[i] + 2.0 * f_[i + 1]) * DX / 6.0;

			// diffusion
		    diag_[i] += 1.0 / DX; 
            right_[i] -= 1.0 / DX;
			diag_[i + 1] += 1.0 / DX;
		    left_[i + 1] -= 1.0 / DX;
        }
	}
}
```

出力ファイル`data_Poisson.txt`からデータを読み取り, gnuplotを用いてグラフを描き, `Poisson.png`に保存する.

```gnuplot
set title '1-dim Poisson equation f=1'

set xl 'x'
set yl 'u'
set grid

plot 'data_Poisson.txt' using 1:2 w l title 'fem 1-dim Poisson eq'

set terminal pngcairo
set output 'Poisson.png'
replot
```

`Poisson.png`は次のようになる. 

![Poisson](https://user-images.githubusercontent.com/75283304/176590279-8a9f36b8-9b57-4c3e-94ec-1b8b19700467.png)

### 1次元定常移流拡散方程式

移流速度$c$で流れ, 拡散係数$D$で拡散する移流拡散方程式は
```math
\frac{\partial u}{\partial t} + \Delta \cdot (cu) 
= \Delta \cdot (D \Delta u)
```
と表せる. 
$1$次元においては
```math
\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} 
= D \frac{\partial^2 u}{\partial t^2}
```
となる. 
定常としたときは
```math
c \frac{du}{dt} = D \frac{d^2 u}{d t^2}
```
となる. 

前回の$1$次元Poisson方程式に移流行列を足し合わせれば良い. 

```cpp
#include <fstream>  // for std::ofstream
#include <iostream> // for std::endl
#include <vector>   // fos std::vector
#include <boost/assert.hpp> // for BOOST_ASSERT

namespace fem{
    class FEM final{
		public:
			using myvec = std::vector<double>;
	
    	private:
        	enum class boundary_condi_type{
            	DIRICLET = 0,
            	LEFT_NEUMANN = 1,
            	RIGHT_NEUMANN = 2
        	};
        
        	static auto constexpr BCT = boundary_condi_type::DIRICLET;

            // coefficient
            static auto constexpr C = 0.1;
            static auto constexpr DIFF = 0.01;

			// boundary conditions
        	static auto constexpr D0 = 0.0; // left Dirichlet
    		static auto constexpr D1 = 1.0; // right Direchlet
			static auto constexpr N0 = 0.5; // left Neumann
	    	static auto constexpr N1 = -1.0; // right Neumann

	    	static auto constexpr ELEMENT = 100;
        	static auto constexpr LENGTH = 1.0;
	    	static auto constexpr NODE = ELEMENT + 1;
    
    	public:
        	static auto constexpr DX = LENGTH / ELEMENT;
	
    	private:
        	myvec bound_;
        	myvec diag_;
			myvec f_;
	    	myvec left_;
	    	myvec right_;

    	public:
        	FEM()
            	: bound_(NODE, 0.0), diag_(NODE, 0.0), f_(NODE), left_(NODE, 0.0), right_(NODE, 0.0)
        	{
				for (auto && elem : f_){
					elem = 0.0; // (right hand side) = 0
				}
			}

        	~FEM() = default;

			// コピーコントラクタでcopy禁止
        	FEM(FEM const &dummy) = delete;

			// operator=()でもcopy禁止
        	FEM & operator=(FEM const &dummy) = delete;

			// boundary condition
        	void boundary();

			// make stiffness matrix
	    	void mat();

			// TDMA method
	    	myvec tdma() const;

			// result output file
			bool output_file(myvec const &x);
    };
}

int main(){
	fem::FEM fem_obj;

	fem_obj.mat();
	fem_obj.boundary();
	auto const x = fem_obj.tdma();

	if (!fem_obj.output_file(x)){
        std::cerr << "output file not open" << std::endl;
        return -1;
    }

	return 0;
}

namespace fem{
    void FEM::boundary(){
		switch (BCT){
            case boundary_condi_type::DIRICLET:
            	diag_[0] = 1.0;
				diag_[NODE - 1] = 1.0;
			    left_[NODE - 1] = 0.0;
        	    right_[0] = 0.0;
				// boundary conditions
				bound_[0] = D0;
        	    bound_[NODE - 1] = D1;
                break;
		
            case boundary_condi_type::LEFT_NEUMANN:
        	    diag_[NODE - 1] = 1.0;
				left_[NODE - 1] = 0.0;
				// boundary conditions
			    bound_[0] -= N0 * DIFF; // left Neumann
				bound_[NODE - 1] = D1;
                break;
			
            case boundary_condi_type::RIGHT_NEUMANN:
			    diag_[0] = 1.0;
        	    right_[0] = 0.0;
				// boundary conditions
        	    bound_[0] = D0;
			    bound_[NODE - 1] += N1 * DIFF; // right Neumann
                break;
            
            default:
                BOOST_ASSERT(!"error! wrong boundary condition");
                break;
        }
	}

	void FEM::mat(){
	    for (auto i = 0; i < ELEMENT; i++){
            // advection
            diag_[i] -= C / 2.0;
            right_[i] += C / 2.0;
            diag_[i + 1] += C / 2.0;
            left_[i + 1] -= C / 2.0; 

			//boundary
		    bound_[i] += (2.0 * f_[i] + 1.0 * f_[i + 1]) * DX / 6.0;
		    bound_[i + 1] += (1.0 * f_[i] + 2.0 * f_[i + 1]) * DX / 6.0;

			// diffusion
		    diag_[i] += DIFF / DX; 
            right_[i] -= DIFF / DX;
			diag_[i + 1] += DIFF / DX;
		    left_[i + 1] -= DIFF / DX; 
        }
	}

	FEM::myvec FEM::tdma() const{
		myvec p(NODE, 0.0);
		myvec q(NODE, 0.0);

		p[0] = -right_[0] / diag_[0];
		q[0] = bound_[0] / diag_[0];

		for (auto i = 1; i < NODE; i++){
		    p[i] = -right_[i] / (diag_[i] + left_[i] * p[i - 1]);
		    q[i] = (bound_[i] - left_[i] * q[i - 1]) / (diag_[i] + left_[i] * p[i - 1]);
	    }

		myvec x(NODE, 0.0);
		x[NODE - 1] = q[NODE - 1];

		for (auto j = NODE - 2; j >= 0; j--){
			x[j] = p[j] * x[j + 1] + q[j];
		}
		
        return x;
 	}

	bool FEM::output_file(myvec const &x){
		std::ofstream ofs("data_advec_diff.txt");
    	if (!ofs) {
        	return false;
    	}

		auto const size = static_cast<int>(x.size());

		for (auto i = 0; i < size; i++){
			ofs << fem::FEM::DX * static_cast<double>(i) << " " << x[i] << std::endl;
		}

		return true;
	}
}
```

出力ファイル`data_advec_diff.txt`からデータを読み込み, gnuplotを用いてグラフを描き, `advec_diff.png`に保存する. 

```gnuplot
set title 'steady-state advection-diffusion equation'

set xl 'x'
set yl 'u'
set grid

plot 'data_advec_diff.txt' using 1:2 w l title 'steady advec diff'

set terminal pngcairo
set output 'steady_advec_diff.png'
replot
```

`advec_diff.png`は次のようになる. 

![steady_advec_diff](https://user-images.githubusercontent.com/75283304/176590809-5b65ff14-f411-4099-a5c8-8520d9f82c5d.png)

### 1次元非定常拡散方程式

拡散係数$D$で拡散する場合の拡散方程式は
```math
\frac{\partial u}{\partial t} = \Delta \cdot (D \Delta u)
```
である. 
$1$次元においては
```math
\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial t^2}
```
となる. 

```cpp
#include <algorithm> // for std::min
#include <array>     // for std::array
#include <fstream>   // for std::ofstream
#include <iostream>  // for std::endl
#include <sstream>   // for std::string, std::stringstream
#include <vector>    // fos std::vector
#include <boost/assert.hpp> // for BOOST_ASSERT
#include <boost/format.hpp> // for boost::format

namespace fem{
    class FEM final{
		public:
			using myvec = std::vector<double>;
	
    	private:
        	enum class boundary_condi_type{
            	DIRICLET = 0,
            	LEFT_NEUMANN = 1,
            	RIGHT_NEUMANN = 2
        	};

            // DIRICLET or LEFT NEUMANN or RIGHT NEUMANN
        	static auto constexpr BCT = boundary_condi_type::RIGHT_NEUMANN;

			// boundary condition
        	static auto constexpr D0 = 0.0; // left Dirichlet
    		static auto constexpr D1 = 0.0; // right Direchlet
			static auto constexpr N0 = 1.0; // left Neumann
	    	static auto constexpr N1 = 0.5; // right Neumann

            // coefficient
            static auto constexpr DIFF = 0.01; // diffusion 
            static auto constexpr THETA = 0.5; // theta method

            // for discrete of x and t
 	    	static auto constexpr ELEMENT = 100; // mesh number
        	static auto constexpr LENGTH = 1.0;  // xrange
            static auto constexpr TEND = 10.0;   // trange
    
    	public:
			static auto constexpr TLOOP = 10000;
            static auto constexpr DT = TEND / TLOOP;
    		static auto constexpr TREP = TLOOP / 1000;
			static auto constexpr NODE = ELEMENT + 1;
			static auto constexpr DX = LENGTH / ELEMENT;
			static auto constexpr KAP = DIFF * DT / (DX * DX);
	
    	private:
        	myvec bound_;
        	myvec diag1_;
			myvec diag2_;
	    	myvec left1_;
			myvec left2_;
	    	myvec right1_;
	    	myvec right2_;
			myvec x_;

    	public:
        	FEM()
            	: bound_(NODE, 0.0), diag1_(NODE, 0.0), diag2_(NODE, 0.0), 
                  left1_(NODE, 0.0), left2_(NODE, 0.0), right1_(NODE, 0.0), right2_(NODE, 0.0), x_(NODE, 0.0)
        	{
				for (auto i = 0; i < NODE; i++){
					x_[i] = std::min(DX * static_cast<double>(i), 1.0 - DX * static_cast<double>(i)); // f = 
				}
			}

        	~FEM() = default;

			// copy constructorでcopy禁止
        	FEM(FEM const &dummy) = delete;

			// operator=()でもcopy禁止
        	FEM & operator=(FEM const &dummy) = delete;

			// boundary condition
        	void boundary();
			void boundary2();

			// TDMA method
	    	void tdma();

			// make stiffness matrix
	    	void mat();

			// output for each TLOOP
			bool result_output(int i);
    };
}

int main(){
    // KAP > 0.5のとき強制終了 
    // Neumannの安定性解析より, KAP <= 0.5でないと正しい結果が得られない
    BOOST_ASSERT(fem::FEM::KAP <= 0.5);

    fem::FEM fem_obj;

	fem_obj.result_output(0);

	if (!fem_obj.result_output(0)){
		std::cerr << "output file not open" << std::endl;
        return -1;
	}

	fem_obj.mat();
	fem_obj.boundary();

	for (auto i = 1; i <= fem::FEM::TLOOP; i++){
		fem_obj.boundary2();
		fem_obj.tdma();

		if (!(i % 1000)){
			if (!fem_obj.result_output(i)){
				std::cerr << "output file not open" << std::endl;
        		return -1;
			}
		}
	}

	return 0;
}

namespace fem{
	bool FEM::result_output(int i){
		std::stringstream ss;
		std::string name;
		std::ofstream ofs;

		ss << i;
		name = ss.str();
		name = "data_nonste_diff_" + name + ".txt";
		
		ofs.open(name.c_str());
		if (!ofs){
			return false;
		}

		for (auto j = 0; j < NODE; j++){
			ofs << DX * static_cast<double>(j) << " " << x_[j] << std::endl;
		}

		return true;
	}

    void FEM::boundary(){
		switch (BCT){
            case boundary_condi_type::DIRICLET:
            	diag1_[0] = 1.0;
				diag1_[NODE - 1] = 1.0;
                diag2_[0] = 1.0;
                diag2_[NODE - 1] = 1.0;
			    left1_[NODE - 1] = 0.0;
                left2_[NODE - 1] = 0.0;
        	    right1_[0] = 0.0;
                right2_[0] = 0.0;
                break;
		
            case boundary_condi_type::LEFT_NEUMANN:
        	    diag1_[NODE - 1] = 1.0;
                diag2_[NODE - 1] = 1.0;
				left1_[NODE - 1] = 0.0;
                left2_[NODE - 1] = 0.0;
                break;
			
            case boundary_condi_type::RIGHT_NEUMANN:
			    diag1_[0] = 1.0;
                diag2_[0] = 1.0;
        	    right1_[0] = 0.0;
                right2_[0] = 0.0;
                break;
            
            default:
                BOOST_ASSERT(!"error! wrong boundary condition");
                break;
        }
	}

    void FEM::boundary2(){
		switch (BCT){
            case boundary_condi_type::DIRICLET:
                bound_[0] = D0;
        	    bound_[NODE - 1] = D1;
                break;
		
            case boundary_condi_type::LEFT_NEUMANN:
                bound_[0] = diag2_[0] * x_[0] + right2_[0] * x_[1] - N0 * DIFF;
				bound_[NODE - 1] = D1;
                break;
			
            case boundary_condi_type::RIGHT_NEUMANN:
                bound_[0] = D0;  
				bound_[NODE - 1] = left2_[NODE - 1] * x_[NODE - 2] + left2_[NODE - 1] * x_[NODE - 1] + N1 * DIFF;
                break;
            
            default:
                BOOST_ASSERT(!"error! wrong boundary condition");
                break;
        }

		for (auto j = 1; j < NODE - 1; j++){
            bound_[j] = left2_[j] * x_[j - 1] + diag2_[j] * x_[j] + right2_[j] * x_[j + 1];
        }
	}

	void FEM::mat(){
	    for (auto i = 0; i < ELEMENT; i++){
            // temporal 1
            diag1_[i] += DX * 2.0 / 6.0 / DT;
            diag1_[i + 1] += DX * 2.0 / 6.0 / DT;
            left1_[i + 1] += DX * 1.0 / 6.0 / DT;
            right1_[i] += DX * 1.0 / 6.0 / DT;

			// diffusion 1
		    diag1_[i] += THETA * DIFF / DX;
            diag1_[i + 1] += THETA * DIFF / DX; 
		    left1_[i + 1] -= THETA * DIFF / DX;
            right1_[i] -= THETA * DIFF / DX;

            // temporal 2
            diag2_[i] += DX * 2.0 / 6.0 / DT;
            diag2_[i + 1] += DX * 2.0 / 6.0 / DT;
            left2_[i + 1] += DX * 1.0 / 6.0 / DT;
            right2_[i] += DX * 1.0 / 6.0 / DT;

			// diffusion 2
		    diag2_[i] -= (1.0 - THETA) * DIFF / DX;
			diag2_[i + 1] -= (1.0 - THETA) * DIFF / DX;
		    left2_[i + 1] += (1.0 - THETA) * DIFF / DX; 
            right2_[i] += (1.0 - THETA) * DIFF / DX;
        }
	}

	void FEM::tdma(){
		myvec p(NODE, 0.0), q(NODE, 0.0);

		p[0] = -right1_[0] / diag1_[0];
		q[0] = bound_[0] / diag1_[0];

		for (auto i = 1; i < NODE; i++){
		    p[i] = -right1_[i] / (diag1_[i] + left1_[i] * p[i - 1]);
		    q[i] = (bound_[i] - left1_[i] * q[i - 1]) / (diag1_[i] + left1_[i] * p[i - 1]);
	    }

		x_[NODE - 1] = q[NODE - 1];

		for (auto i = NODE - 2; i >= 0; i--){
			x_[i] = p[i] * x_[i + 1] + q[i];
		}
 	}
}
```

出力ファイル`data_nonste_diff.txt`からデータを読み込み, gnuplotを用いてグラフを描き, `nonste_diff.png`に保存する. 

```gnuplot
set title 'non-steady diffusion eq'

set xl 'x'
set yl 'u'
set grid

plot 'data_nonste_diff_0.txt'       using 1:2 w l title 't = 0'
replot 'data_nonste_diff_1000.txt'  using 1:2 w l title 't = 1'
replot 'data_nonste_diff_3000.txt'  using 1:2 w l title 't = 3'
replot 'data_nonste_diff_5000.txt'  using 1:2 w l title 't = 5'
replot 'data_nonste_diff_7000.txt'  using 1:2 w l title 't = 7'
replot 'data_nonste_diff_9000.txt'  using 1:2 w l title 't = 9'
replot 'data_nonste_diff_10000.txt' using 1:2 w l title 't = 10'

set terminal pngcairo
set output 'nonste_diff.png'
replot
```

`nonste_diff.png`は次のようになる. 

![nonste_diff](https://user-images.githubusercontent.com/75283304/176592047-7c408c39-5438-4d39-a9df-3f29f7b33abb.png)

### 1次元非定常移流拡散方程式

移流速度$c$で流れ, 拡散係数$D$で拡散する場合の移流拡散方程式は
```math
\frac{\partial u}{\partial t} + \Delta \cdot (cu) 
= \Delta \cdot (D \Delta u)
```
である. 
$1$次元においては
```math
\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} 
= D \frac{\partial^2 u}{\partial t^2}
```
となる. 

```cpp
#include <algorithm> // for std::min
#include <array>     // for std::array
#include <fstream>   // for std::ofstream
#include <iostream>  // for std::endl
#include <sstream>   // for std::string, std::stringstream
#include <vector>    // fos std::vector
#include <boost/assert.hpp> // for BOOST_ASSERT
#include <boost/format.hpp> // for boost::format

namespace fem{
    class FEM final{
		public:
			using myvec = std::vector<double>;
	
    	private:
        	enum class boundary_condi_type{
            	DIRICLET = 0,
            	LEFT_NEUMANN = 1,
            	RIGHT_NEUMANN = 2
        	};

            // DIRICLET or LEFT NEUMANN or RIGHT NEUMANN
        	static auto constexpr BCT = boundary_condi_type::DIRICLET;

			// boundary condition
        	static auto constexpr D0 = 0.0; // left Dirichlet
    		static auto constexpr D1 = 0.0; // right Direchlet
			static auto constexpr N0 = 1.0; // left Neumann
	    	static auto constexpr N1 = 0.5; // right Neumann

            // coefficient
			static auto constexpr C = 0.1;     // advection
            static auto constexpr DIFF = 0.01; // diffusion 
            static auto constexpr THETA = 0.5; // theta method

            // for discrete of x and t
 	    	static auto constexpr ELEMENT = 100; // mesh number
        	static auto constexpr LENGTH = 1.0;  // xrange
            static auto constexpr TEND = 10.0;   // trange
    
    	public:
			static auto constexpr TLOOP = 10000;
            static auto constexpr DT = TEND / TLOOP;
    		static auto constexpr TREP = TLOOP / 1000;
			static auto constexpr NODE = ELEMENT + 1;
			static auto constexpr DX = LENGTH / ELEMENT;
			static auto constexpr KAP = DIFF * DT / (DX * DX);
	
    	private:
        	myvec bound_;
        	myvec diag1_;
			myvec diag2_;
	    	myvec left1_;
			myvec left2_;
	    	myvec right1_;
	    	myvec right2_;
			myvec x_;

    	public:
        	FEM()
            	: bound_(NODE, 0.0), diag1_(NODE, 0.0), diag2_(NODE, 0.0), 
                  left1_(NODE, 0.0), left2_(NODE, 0.0), right1_(NODE, 0.0), right2_(NODE, 0.0), x_(NODE, 0.0)
        	{
				for (auto i = 0; i < NODE; i++){
					x_[i] = std::min(DX * static_cast<double>(i), 1.0 - DX * static_cast<double>(i)); // f = 
				}
			}

        	~FEM() = default;

			// copy constructorでcopy禁止
        	FEM(FEM const &dummy) = delete;

			// operator=()でもcopy禁止
        	FEM & operator=(FEM const &dummy) = delete;

			// output for each TLOOP
			bool result_output(int i);

			// boundary condition
        	void boundary();
			void boundary2();

			// TDMA method
	    	void tdma();

			// make stiffness matrix
	    	void mat();
    };
}

int main(){
    // KAP > 0.5のとき強制終了 
    // Neumannの安定性解析より, KAP <= 0.5でないと正しい結果が得られない
    BOOST_ASSERT(fem::FEM::KAP <= 0.5);

    fem::FEM fem_obj;

	fem_obj.result_output(0);

	if (!fem_obj.result_output(0)){
		std::cerr << "output file not open" << std::endl;
        return -1;
	}

	fem_obj.mat();
	fem_obj.boundary();

	for (auto i = 1; i <= fem::FEM::TLOOP; i++){
		fem_obj.boundary2();
		fem_obj.tdma();

		if (!(i % 1000)){
			if (!fem_obj.result_output(i)){
				std::cerr << "output file not open" << std::endl;
        		return -1;
			}
		}
	}

	return 0;
}

namespace fem{
	bool FEM::result_output(int i){
		std::stringstream ss;
		std::string name;
		std::ofstream ofs;

		ss << i;
		name = ss.str();
		name = "data_nonste_advec_diff_" + name + ".txt";
		
		ofs.open(name.c_str());
		if (!ofs){
			return false;
		}

		for (auto j = 0; j < NODE; j++){
			ofs << DX * static_cast<double>(j) << " " << x_[j] << std::endl;
		}

		return true;
	}

    void FEM::boundary(){
		switch (BCT){
            case boundary_condi_type::DIRICLET:
            	diag1_[0] = 1.0;
				diag1_[NODE - 1] = 1.0;
                diag2_[0] = 1.0;
                diag2_[NODE - 1] = 1.0;
			    left1_[NODE - 1] = 0.0;
                left2_[NODE - 1] = 0.0;
        	    right1_[0] = 0.0;
                right2_[0] = 0.0;
                break;
		
            case boundary_condi_type::LEFT_NEUMANN:
        	    diag1_[NODE - 1] = 1.0;
                diag2_[NODE - 1] = 1.0;
				left1_[NODE - 1] = 0.0;
                left2_[NODE - 1] = 0.0;
                break;
			
            case boundary_condi_type::RIGHT_NEUMANN:
			    diag1_[0] = 1.0;
                diag2_[0] = 1.0;
        	    right1_[0] = 0.0;
                right2_[0] = 0.0;
                break;
            
            default:
                BOOST_ASSERT(!"error! wrong boundary condition");
                break;
        }
	}

    void FEM::boundary2(){
		switch (BCT){
            case boundary_condi_type::DIRICLET:
                bound_[0] = D0;
        	    bound_[NODE - 1] = D1;
                break;
		
            case boundary_condi_type::LEFT_NEUMANN:
                bound_[0] = diag2_[0] * x_[0] + right2_[0] * x_[1] - N0 * DIFF;
				bound_[NODE - 1] = D1;
                break;
			
            case boundary_condi_type::RIGHT_NEUMANN:
                bound_[0] = D0;  
				bound_[NODE - 1] = left2_[NODE - 1] * x_[NODE - 2] + left2_[NODE - 1] * x_[NODE - 1] + N1 * DIFF;
                break;
            
            default:
                BOOST_ASSERT(!"error! wrong boundary condition");
                break;
        }

		for (auto j = 1; j < NODE - 1; j++){
            bound_[j] = left2_[j] * x_[j - 1] + diag2_[j] * x_[j] + right2_[j] * x_[j + 1];
        }
	}

	void FEM::mat(){
	    for (auto i = 0; i < ELEMENT; i++){
            // temporal 1
            diag1_[i] += DX * 2.0 / 6.0 / DT;
            diag1_[i + 1] += DX * 2.0 / 6.0 / DT;
            left1_[i + 1] += DX * 1.0 / 6.0 / DT;
            right1_[i] += DX * 1.0 / 6.0 / DT;

			//advection 1
			diag1_[i] -= THETA * C / 2.0;
			diag1_[i + 1] += THETA * C / 2.0;
			left1_[i + 1] -= THETA * C / 2.0;
			right1_[i] += THETA * C / 2.0;

			// diffusion 1
		    diag1_[i] += THETA * DIFF / DX;
            diag1_[i + 1] += THETA * DIFF / DX; 
		    left1_[i + 1] -= THETA * DIFF / DX;
            right1_[i] -= THETA * DIFF / DX;

            // temporal 2
            diag2_[i] += DX * 2.0 / 6.0 / DT;
            diag2_[i + 1] += DX * 2.0 / 6.0 / DT;
            left2_[i + 1] += DX * 1.0 / 6.0 / DT;
            right2_[i] += DX * 1.0 / 6.0 / DT;

			//advection 2
			diag2_[i] += (1 - THETA) * C / 2.0;
			diag2_[i + 1] -= (1 - THETA) * C / 2.0;
			left2_[i + 1] += (1 - THETA) * C / 2.0;
			right2_[i] -= (1 - THETA) * C / 2.0;

			// diffusion 2
		    diag2_[i] -= (1.0 - THETA) * DIFF / DX;
			diag2_[i + 1] -= (1.0 - THETA) * DIFF / DX;
		    left2_[i + 1] += (1.0 - THETA) * DIFF / DX; 
            right2_[i] += (1.0 - THETA) * DIFF / DX;
        }
	}

	void FEM::tdma(){
		myvec p(NODE, 0.0), q(NODE, 0.0);

		p[0] = -right1_[0] / diag1_[0];
		q[0] = bound_[0] / diag1_[0];

		for (auto i = 1; i < NODE; i++){
		    p[i] = -right1_[i] / (diag1_[i] + left1_[i] * p[i - 1]);
		    q[i] = (bound_[i] - left1_[i] * q[i - 1]) / (diag1_[i] + left1_[i] * p[i - 1]);
	    }

		x_[NODE - 1] = q[NODE - 1];

		for (auto i = NODE - 2; i >= 0; i--){
			x_[i] = p[i] * x_[i + 1] + q[i];
		}
 	}
}
```

出力ファイル`data_nonste_advec_diff.txt`からデータを読み込み, gnuplotを用いてグラフを描き, `nonste_advec_diff.png`に保存する. 

```gnuplot
set title 'non-steady advection diffusion eq'

set xl 'x'
set yl 'u'
set grid

plot 'data_nonste_advec_diff_0.txt'       using 1:2 w l title 't = 0'
replot 'data_nonste_advec_diff_1000.txt'  using 1:2 w l title 't = 1'
replot 'data_nonste_advec_diff_3000.txt'  using 1:2 w l title 't = 3'
replot 'data_nonste_advec_diff_5000.txt'  using 1:2 w l title 't = 5'
replot 'data_nonste_advec_diff_8000.txt'  using 1:2 w l title 't = 7'
replot 'data_nonste_advec_diff_10000.txt' using 1:2 w l title 't = 10'

set terminal pngcairo
set output 'nonste_advec_diff.png'
replot
```

`nonste_advec_diff.png`は次のようになる.

![nonste_advec_diff](https://user-images.githubusercontent.com/75283304/176592661-7c3b45ca-2d2e-4239-9b32-b9c6b38e56f1.png)

### 1次元Burgers方程式

1次元のBurgers方程式は
```math
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} 
= D \frac{\partial^2 u}{\partial t^2}
```
となる. 
左辺第2項目が非線形項であるが, Picard反復を用いて線形化する. 
理論的にはHopf-Cole変換によって拡散方程式に帰着させることができる. 

## その他の解法

(有限)差分法や有限要素法のほかにも有限体積法がある. 

# 参考文献

全体を通して

+ 小高知宏 「Cによる数値計算とシミュレーション」
+ [数値解析についての資料](http://kobayashi.hub.hit-u.ac.jp/lecture/na.pdf)

連立一次方程式に関して

+ [mk-mode BLOG](https://komasaru.github.io/)
+ [数学とか語学とか楽しいよね](https://mathlang.hatenablog.com/)

常微分方程式に関して

+ [微分方程式の数値解法](http://www.den.t.u-tokyo.ac.jp/ad_prog/ode/)
+ [常微分方程式の初期値問題を解くプログラムの書き方](http://nalab.mind.meiji.ac.jp/~mk/labo/text/intro-ode-simulation.pdf)

偏微分方程式に関して

+ [有限要素法プログラミング演習](https://web.archive.org/web/20180922014124/http://www.sml.k.u-tokyo.ac.jp/members/nabe/FEM/FEM.pdf)
+ [工学系のための偏微分方程式 数値解法のＣ言語による実習](https://www.sci.kumamoto-u.ac.jp/physics/plasma/jisshusho.pdf)
+ [基礎数値解析 -偏微分方程式の数値シミュレーション技法入門-](https://www1.gifu-u.ac.jp/~tanaka/numerical_analysis.pdf)
+ [数学とか語学とか楽しいよね](https://mathlang.hatenablog.com/)

Eigenに関して

+ [C++行列計算ライブラリEigen入門](https://qiita.com/yohm/items/a03006790dc1e54a87be)
+ [Eigen ー C++で線形代数を！(3)](http://www.singularpoint.org/blog/math/eigen-3/)
+ [Eigen - C++で使える線形代数ライブラリ](http://blog.livedoor.jp/tek_nishi/archives/8623876.html)

boostに関して

+ [Numerical integration in C++ with Boost odeint](http://boccelliengineering.altervista.org/junk/boost_integration/boost_odeint.html)
+ [C++で常微分方程式：boost::odeint入門](https://qiita.com/hmito/items/483445ac0d42fb4428a5)
+ [Boost.odeintの紹介](https://qiita.com/ignis_fatuus/items/74cb32815753d891a716#:~:text=Boost.Odeint%E3%82%92%E4%BD%BF%E3%81%A3%E3%81%A6%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%82%92%E8%A7%A3%E3%81%8F%E9%83%A8%E5%88%86%E3%82%92%E6%9B%B8%E3%81%8D%E8%B6%B3%E3%81%97%E3%81%BE%E3%81%99%E3%80%82%20Integrate%E3%81%AF%E7%A9%8D%E5%88%86%E3%81%99%E3%82%8B%E3%81%A8%E3%81%84%E3%81%86%E6%84%8F%E5%91%B3%E3%81%A7%E3%81%99%E3%81%8C%E3%80%81%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%82%92%E8%A7%A3%E3%81%8F%20%28x%20%E3%82%92%E6%B1%82%E3%82%81%E3%82%8B%29%E3%81%A8%E3%81%84%E3%81%86%E3%81%AE%E3%81%AF%E7%B5%90%E5%B1%80%20x%20%3D%20%E2%88%AB,d%20x%20d%20t%20d%20t%20%E3%82%92%E8%A8%88%E7%AE%97%E3%81%99%E3%82%8B%E3%81%93%E3%81%A8%E3%81%AA%E3%81%AE%E3%81%A7%E3%80%81%E7%A9%8D%E5%88%86%E3%81%99%E3%82%8B%E3%81%93%E3%81%A8%E3%81%A8%E5%90%8C%E3%81%98%E3%81%A7%E3%81%99%E3%80%82)

gnuplotに関して
+ [gnuplot の基本的な使い方](https://qiita.com/python_walker/items/2b4febdba01c2a679169)