数値計算の基本をC++で書いていきます.
YonohaのGithubにも全体のコードとグラフを公開しています.
boostやEigen, gnuplotを使用しているコードのコンパイルには, boostやEigen, gnuplotのインストールがそれぞれ必要です.
変更した方が良い箇所があれば, この記事かTwitterの@Yonono_01まで連絡ください.
目次
偏微分方程式 (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に関して