数値計算の基本をC++で書いていきます.
YonohaのGithubにも全体のコードとグラフを公開しています.

boostやEigen, gnuplotを使用しているコードのコンパイルには, boostやEigen, gnuplotのインストールがそれぞれ必要です.

変更した方が良い箇所があれば, この記事かTwitterの@Yonono_01まで連絡ください.

目次

偏微分方程式 (Partial Differential Equation)
参考文献

偏微分方程式 (Partial Differential Equation)

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

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

差分法 (Difference Method)

次の$1$次元拡散方程式

\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} ~ (0 < x < 1, 0 < t < \infty)

を初期条件

u(x,0) = \sin(\pi x) ~ (0 \leq x \leq 1)

と境界条件

u(0,t) = u(1,t) = 0 ~ (0 < t < \infty)

の元で解く.
この偏微分方程式をメッシュ数を$I$として, 差分スキームの形に書きなおすと

差分方程式

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)

初期条件

u_i^0 = \sin(\pi x_i)

境界条件

u_0^n = u_I^n = 0

のようになる.

#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に保存する.

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は次のようになる.

有限要素法 (Finite Element Method)

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

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

1次元Poisson方程式

$f$を既知関数とし, $u$を求める未知関数とする.
このとき, 次の偏微分方程式

\Delta u + f = 0

をPoisson方程式という.
$1$次元においては

\frac{d^2 u}{dx^2} + f = 0 ~ (x \in [0, L])

となる.
今回は

f = 1

の場合を考える.

#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に保存する.

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は次のようになる.

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

移流速度$c$で流れ, 拡散係数$D$で拡散する移流拡散方程式は

\frac{\partial u}{\partial t} + \Delta \cdot (cu) 
= \Delta \cdot (D \Delta u)

と表せる.
$1$次元においては

\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} 
= D \frac{\partial^2 u}{\partial t^2}

となる.
定常としたときは

c \frac{du}{dt} = D \frac{d^2 u}{d t^2}

となる.

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

#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に保存する.

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は次のようになる.

1次元非定常拡散方程式

拡散係数$D$で拡散する場合の拡散方程式は

\frac{\partial u}{\partial t} = \Delta \cdot (D \Delta u)

である.
$1$次元においては

\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial t^2}

となる.

#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に保存する.

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は次のようになる.

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

移流速度$c$で流れ, 拡散係数$D$で拡散する場合の移流拡散方程式は

\frac{\partial u}{\partial t} + \Delta \cdot (cu) 
= \Delta \cdot (D \Delta u)

である.
$1$次元においては

\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} 
= D \frac{\partial^2 u}{\partial t^2}

となる.

#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に保存する.

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は次のようになる.

1次元Burgers方程式

1次元のBurgers方程式は

\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} 
= D \frac{\partial^2 u}{\partial t^2}

となる.
左辺第2項目が非線形項であるが, Picard反復を用いて線形化する.
理論的にはHopf-Cole変換によって拡散方程式に帰着させることができる.

その他の解法

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

参考文献

全体を通して

小高知宏「Cによる数値計算とシミュレーション」
数値解析についての資料

連立一次方程式に関して

常微分方程式に関して

偏微分方程式に関して

Eigenに関して

boostに関して

gnuplotに関して

gnuplot の基本的な使い方