偏微分方程式I

こんな感じ：

// main.cpp
#include <iostream>
#include "CONV.hpp"
int main()
{
    size_t nx=100; // 空間分割数
    double dt=0.01; // 時間刻み幅
    size_t ntmax=1000; // 最大時間ステップ数
    size_t ntsave=10; // データ保存間隔
    CONV mySolver(nx,dt,"Data/pde100.data");
    mySolver.Write();                  //初期値も保存しよう
    for(size_t nt=0;nt<ntmax;nt++)
    {
        mySolver.Step();
        if (nt % ntsave == 0) mySolver.Write();
    }
    return 0;
}

定義の読み込みと, メンバー変数のところは：

//  CONV.hpp

#pragma once
#include <iostream>
#include <filesystem>
#include <cmath>
#include <vector>
#include <fstream>
class CONV
{
   double dt, CFL;
   double time;
   class Point  //めんどくさかったので, 内部クラスにした
   {
   public:
       size_t id;                // これは・・・いらないかもな？
       double X;
       double U_VAL; // U_n
       double U_PRE; // U_{n-1}
   }; 
   std::vector<Point> Data;
   std::ofstream ofs;

コンストラクタは, 保存するファイル名を入れて初期化する予定でしたが．なーんか，それ本質じゃ無い気がしてきて

• ほれ．コンストラクタは「最初にやること」ですよね．できたファイルを保存する名前，は，そうじゃ無い気が
• むしろ「初期値を指定する」のが良いのでは無いか？とおもったので, 仕様を変更．ラムダ式で与えても良いが，忘れるユーザの対策もしておいた

public:
   CONV(const size_t& nx, const double& dt, const std::function<double(double)>& U0 = [](double x)
                                          { return 0.; }) : dt(dt)
   {
       double xmax = 1.0;
       time = 0.0;            // 初期時間
       double dx = xmax / nx; // 空間刻み幅の計算
       CFL = dt / dx;         // CFL条件の計算
       Data.resize(nx + 1);   // 空間分割数に基づいて初期化
       for (size_t i = 0; i < Data.size(); ++i) Data[i].id = i; // 各点のIDを設定
       for (auto& P:Data) P.U_PRE = P.U_VAL = U0(P.X = P.id * dx); // 初期条件の適用
       std::cout << "Starting CONV solver with nx=" << Data.size() << ", dt=" << dt << std::endl;
    }

おや．ここいらへんの const とは，なんでしょう？　それは後で説明しますね！

書き出しファイル名は, Write()を最初に呼び出すときに与え, それ以外では引数なしにしましょう．

bool Write(const std::string& filename)
   {
       if (ofs) ofs.close(); // 既存のファイルストリームを閉じる
       std::filesystem::path datafile = getenv("HOME");
       datafile /= filename;　// 引数のファイル名は，$HOME に相対的なモニに限定
       std::filesystem::create_directories(datafile.parent_path());
       ofs.open(datafile);
       if (!ofs)　throw std::ios_base::failure("Failed to open file: " + datafile.string());
       return Write();       // 引数なし版: 現在のofsにデータを書き出す
   }
   bool Write()
   const {
       try
       {
           ofs << std::endl << "# time= " << time << std::endl;
           for (auto &v : Data) ofs << v.X << " " << v.U_VAL << std::endl;
       }
       catch (const std::exception &e)
       {
           std::cerr << "Error writing to file: " << e.what() << std::endl;
           return false;
       }
       return true;
   }

Stepはスキーム通りです. Step()の引数に，境界条件の該当時刻の値を与える，という仕様訂正をしました．

   void Step(const double U_W=1.0)
   {
       time+=dt;
       Data[0].U_VAL = U_W;
       for(size_t i=1;i<Data.size();++i)
           Data[i].U_VAL=(1.-CFL)*Data[i].U_PRE + CFL*Data[i-1].U_PRE;
       for(auto& v:Data) v.U_PRE=v.U_VAL;
   }

   double get_time() const { return time; }

数値計算的には, この辺で, 計算した分布を, 以前の分布として交換するのを忘れるとまずいですよね！この関数の意味は，後で説明しますね！

void Step_implicitと書いたら, AIモンキーが書きましたが，間違ってたので直させました. ううん，結構AI嘘つきますね・・・2030年くらいまでは用心して使いましょう.

   void Step_Implicit(const double U_W=1.0)
   {
       time+=dt;
       Data[0].U_VAL = U_W;
       for(size_t i=1;i<Data.size();++i) Data[i].U_VAL = (Data[i].U_PRE+CFL*Data[i-1].U_VAL)/(1.+CFL);
       for(auto& v:Data) v.U_PRE=v.U_VAL;
   }

mainの方は

// main.cpp
#include <iostream>
#include "CONV.hpp"
int main()
{
   size_t nx=10; // 空間分割数
   double dt=0.1; // 時間刻み幅
   size_t ntmax=50; // 最大時間ステップ数
   size_t ntsave=10; // データ保存間隔
   CONV mySolver(nx,dt);
   mySolver.Write("Data/test03/PDE_02.dat"); // 初期状態の保存
   for(size_t nt=0;nt<ntmax;nt++)
   {
       // mySolver.Step();       // エクスプリシットでの時間ステップ
       mySolver.Step_Implicit(); // インプリシットでの時間ステップ
       if ((nt+1) % ntsave == 0) mySolver.Write();
   }
   return 0;
}

結果を比較してみると:
やっぱさ，なんか，インプリシットって，ダレるんだよな

#!/bin/bash
#   PDE1.sh: 偏微分方程式I
. draw.sh
. standards.sh
COLOR=1
X_AXIS='\fIX\sb\fR1\eb'
X_RANGE='0 10 2.5'
Y_AXIS='\fIu'
Y_RANGE='-0.25 1.25 0.25'
LINE=( Solid:Black:0.005 Solid:Red:0.003 )
SYMBOL=( Omit )
cat > $TMP/legend.dat <<@
#analytical solution t=1
#m=Dash:Red,S=Omit     ←PlotUtils流力研PatchLevel9では，凡例で線の太さが変えられない
0 1
1 1
1 0
10 0

#analytical solution t=2
#m=Dot:Blue,S=Omit
0 1
2 1
2 0
10 0

#m=Solid:Black,S=Omit
3.25 0.5
5.0 1.0  lb:\fIt\fR = 3

#m=Solid:Red,S=Omit
6 0.65
7 0.65  l: implicit

#m=Solid:Black,S=Omit
6 0.75
7 0.75  l: explicit
@
LEGEND=$TMP/legend.dat
draw my.pdf PDE_01.dat PDE_02.dat

うえは draw.sh でグラフを描いた例です．こんな簡単な図の修正に300秒以上必要なソフトは不良品ですので，消費者センターに通報しましょう

できたつもり（バグっている）

// pvector.hpp
#pragma once
#include <iostream>
#include <vector>
#include <cmath>
class pvector
{
   std::vector<double> data;
public:
   pvector(size_t size = 0) : data(size) {}
   auto &operator[](size_t index)
   {return data[index % data.size()];}
   const auto &operator[](size_t index) const
   {return data[index % data.size()];}
   auto size() const
   {return data.size();}
   void resize(size_t newSize)
   {data.resize(newSize);}
   // range-based for 用イテレータ
   auto begin() { return data.begin(); }
   auto end() { return data.end(); }
   auto begin() const { return data.begin(); }
   auto end() const { return data.end(); }
};

お客様に届ける前に，動作点検しましょうね．

// main.cpp
#include <iostream>
...
#define _TEST
#ifndef _TEST
...
#else
   std::cout << "This is a test build. No PDE solver executed." << std::endl;
   pvector Data(7);
   std::cout << "SIZE: " << Data.size() << std::endl;
   for (size_t n=0;auto &v : Data) v=n++;
   for(int i=-10;i<10;i++) std::cout << "Data[" << i << "] = " << Data[i] << std::endl;
#endif
...

ここは範囲forで, 数えることができるという小技です．まあどうでも良いか．いや, ここで, n++ と ++n は異なる．というのもポイントですかね．

• 前者はnをvに代入してから, nを1増やす．Data[0]=0, Data[1]=1, ... になる

• 後者は, nを1増やしてから, vに代入する．Data[0]=1. Data[2]=2.... になる

問題は，ここで，Data[-3]だとどうなるんだ？ということです．クラスの方では引数がsize_tですので，負の数値を受け取れません．ここは正しくはintにしとかないと，ダメです．それで良いのかと思って実行すると，あらあら不思議，

This is a test build. No PDE solver executed.
SIZE: 7
...
Data[-3] = 6　　　←ありゃ？ここは，4
Data[-2] = 0　　　←ありゃ？ここは，5
Data[-1] = 1　　　←ありゃ？ここは，6
Data[0] = 0
Data[1] = 1
...

C++のN%Mの実装は, Mが整数の時と非負整数の時で異なり，ちょっと訳がわからないです．

       int n1 = 7;
       unsigned int n2 = 7;
       size_t n3= 7;              // index=-1  -2  -3  -4 ... -6  -7
       int val1 = index % n1;     // val1 =-1  -2  -3  -4 ... -6   0  これは   sign(N)*abs(N)%M
       int val2 = index % n2;     // val2 = 3   2   1   0      5   3  おそらく -7 を32bit非負整数に変換してから計算
       int val3 = index % n3;     // val3 = 1   0   6   5      3   2　おそらく -7 を64bit非負整数に変換してから計算

こんなものと格闘している暇はないので，AIモンキーに直してもらいました. あれ？僕が使ってるのより式が簡単だな

// pvector.hpp
#pragma once
#include <iostream>
#include <vector>
#include <cmath>
class pvector
{
   std::vector<double> data;
   auto pindex(int index) // インデックスを正の値に変換するヘルパー関数
   const {
       int n = data.size();  // size()は非負整数size_tなのでintに変更
       return ((index % n) + n) % n;
   }
public:
   pvector(size_t size = 0) : data(size) {}    
   auto &operator[](int index)
   {return data[pindex(index)];}
   const auto &operator[](int index) const
   {return data[pindex(index)];}
   auto size() const
   {return data.size();}
   void resize(size_t newSize)
   {data.resize(newSize);}
   auto begin() { return data.begin(); }
   auto end() { return data.end(); }
   auto begin() const { return data.begin(); }
   auto end() const { return data.end(); }
   auto front() { return data.front(); }
   auto back() { return data.back(); }
};

この辺は，後の都合で入れてあります．これで正常: 

This is a test build. No PDE solver executed.
SIZE: 7
Data[-10] = 4
Data[-9] = 5
Data[-8] = 6
Data[-7] = 0
Data[-6] = 1
Data[-5] = 2
Data[-4] = 3
Data[-3] = 4
Data[-2] = 5
Data[-1] = 6
Data[0] = 0
Data[1] = 1
Data[2] = 2
Data[3] = 3
Data[4] = 4
Data[5] = 5
Data[6] = 6
Data[7] = 0
Data[8] = 1
Data[9] = 2

♪ 

下がtemplate です．template はTを与えるとクラスになるものであって，Tが未定の場合，クラスではないです．

#pragma once
#include <iostream>
#include <vector>
#include <cmath>
template<typename T>
class pvector
{
   std::vector<T> data;
   auto pindex(int index) // インデックスを正の値に変換するヘルパー関数
   const {
       int n = data.size();
       return ((index % n) + n) % n;
   }
public:
   pvector(size_t size = 0) : data(size) {}    
   auto &operator[](int index) {return data[pindex(index)];}
   const auto &operator[](int index) const {return data[pindex(index)];}
   auto size() const {return data.size();}
   void resize(size_t newSize) {data.resize(newSize);}
   auto begin() { return data.begin(); }
   auto end() { return data.end(); }
   auto begin() const { return data.begin(); }
   auto end() const { return data.end(); }
};

変更されたのは黄色のところだけです．Tのところがdoubleならdoubleの配列が，myClassならmyClassの配列ができます．もちろん，不用意に T と書いてある型名は，全て検索置換されてしまうので，気をつけましょうね．

まずは概要ですよね.

// pCONV.hpp    Solve advection equation using explicit method. Periodic problem
#pragma once
#include <iostream>
#include <filesystem>
#include <cmath>
#include <vector>
#include <fstream>
#include "pvector.hpp"
class pCONV
{
   double dt, CFL;
   double time;
   int nt = 0; // 時間ステップのカウンタ
   class Point
   {
   public:
       Point* pre = nullptr;   // ポインタを追加して前のポイントへのリンクを作成
       double X;
       pvector<double> U;      // U_n
       Point() : U(3) {}       // デフォルトコンストラクタでサイズ3のpvectorを初期化
       auto &operator[](int index) { return U[index]; }
       const auto &operator[](int index) const { return U[index]; }
   };
   pvector<Point> Data;
   std::ofstream ofs;
public:
   pCONV(const size_t &nx, const double &dt, const std::function<double(double)> &U0 = [](double x)
                                            { return 0.; }) : dt(dt)
   {
       double xmax = 1.0;
       time = 0.0;             // 初期時間
       double dx = xmax / nx;  // 空間刻み幅の計算
       CFL = dt / dx;          // CFL条件の計算
       Data.resize(nx);        // 空間分割数に基づいて初期化
       for (size_t i = 0; auto &P : Data)
       {
           P.pre = &Data[i - 1];
           P[nt] = P[nt - 1] = P[nt - 2] = U0(P.X = i * dx); // 初期条件の適用
           i++;
       }
       std::cout << "Starting pCONV solver with nx=" << Data.size() << ", dt=" << dt << std::endl;
       for(auto& P:Data)
       {
           std::cout << std::setprecision(5) << std::fixed << std::showpos 
           << P.X << " : " << P[nt] << " " << (*P.pre)[nt] << std::endl;
       }
   }

ここは新顔です．前回, Pointに番号をつけたのですが, あまり使い道がなかったので，今回は左側のPointへのポインターを装備してみました．Pointクラスは，ここではまだ定義中であるのですが，自分自身へのポインタは使えるのです．

そもそも，ポインターは「住所」ですので，「あー住所か．なら全角文字30文字の入力欄作っとけば，まあ，はいるやろ」という程度しか考える必要がないので，

class B;
class Point {
...
    B* ptrB;
...

と書いても問題はない．Bは少なくともクラスであることは明らかで，ここには住所欄を作れば良いだけ．「Bに含まれるこの関数」「Bに含まれるこの変数」とかを述べるのでなければ，OKです．こういうのを, 前方宣言と言います．

で，prevフィールドの使い方ですが，ここで，左隣の点を指差せ，と指定していますね！

Pointクラスにもオペレータを装備しているので, Data[nx]が, nx番目の格子点のデータで, 時刻ntのUの値は, Data[nx][nt] でアクセスできることになりました．で，この辺はテストです. 

for(auto& P:Data)で, PはPointインスタンスになっている. P.preは左隣のPointを指さしているので＊でインスタンスに戻して, それも値を画面出力．このへんは, 桁数を固定とか, 符号出しとけ，という指示で，

Starting pCONV solver with nx=10, dt=0.02
+0.00000 : +0.00000 -0.07016
+0.10000 : +0.07016 +0.00000
+0.20000 : +0.77809 +0.07016
+0.30000 : +0.77809 +0.77809
+0.40000 : +0.07016 +0.77809
+0.50000 : +0.00000 +0.07016
+0.60000 : -0.07016 +0.00000
+0.70000 : -0.77809 -0.07016
+0.80000 : -0.77809 -0.77809
+0.90000 : -0.07016 -0.77809

ほれ．Xの値 : 自分の値と，左隣の値が出力されてますね！これは, main()で

pCONV mySolver(nx,dt,[](double x) { return std::pow(std::sin(2.0 * M_PI * x),5); }); // 初期条件として正弦波を設定

と指定したケースです．初期値が\(\sin^5(2\pi x)\)となっています.

では時間発展を記述しましょう. それ以外のWriteとかは，流用で．前回作成した1次精度・妖怪法スキームは

//  pCONV.hpp

...

   void StepLQ()   //  1次精度, explicit method

   {

       time = (++nt) * dt;

       for(auto& P:Data) P[nt]=(1.-CFL)*P[nt-1]+ CFL*(*P.pre)[nt-1];

   }

こうなっちゃいますね. Pは各格子点データ, P[nt]は(ここで1増やしているので)次時刻の値で，これは左の格子点のデータです. では, 空間微分を高次精度にしてみましょう:

   void StepHQ()
   {
       time = (++nt) * dt;
       for (auto &P : Data)
       {
           auto& Pn1=*P.pre;   //左隣
           auto& Pn2=*Pn1.pre; //左左隣
           P[nt] = (1.-1.5*CFL)*P[nt - 1] +2.*CFL*Pn1[nt-1]-0.5*CFL* Pn2[nt - 1];
       }
   }

結果はこんな感じで　緑が厳密解です．格子点数10, CFL=0.2

低次精度では1周する間に波が減衰してしまっている．時間1次-空間2次のは，波の減衰は少ないが，波の速度が，あまり正確ではないですね・・・

• 空間格子数が10点では，流石に難しいな．もう少し増やしてみよう．
• 他のスキームも試してみたい・・・となると，次の継承を覚えた方が良い．
• なお，この図を描くdrawスクリプトは

#!/bin/bash
#   PDEh.sh: 偏微分方程式I: 周期解
. draw.sh
. standards.sh
COLOR=1
X_AXIS='\fIX\sb\fR1\eb'
X_RANGE='0 1 0.25'
Y_AXIS='\fIu'
Y_RANGE='-1.25 1.25 0.5'
LINE=( Solid:Black:0.002 Solid:Red:0.003 Solid:Green:0.003 )
SYMBOL=( Omit )
#解析解のデータが，欲しい
seq 0 1|\
   gawk '
BEGIN{pi=4*atan2(1,1);}
{
time=$1*0.5;
printf("#time= %f\n",time);
for(i=0;i<101;i++)
{
   x=i/100;y=sin(2*pi*(x+time))^5;printf("%f %f\n",x,y);
}
printf("\n");
}
END{printf("\n");}' > $TMP/sol

#凡例
cat > $TMP/legend.dat <<@
#m=Solid:Black,S=Omit
0.1 1.1
0.2 1.1 l:low prec.
#m=Solid:Red,S=Omit
0.6 1.1
0.7 1.1 l:high prec.
#m=Solid:Black,S=Omit
0.70 0.25 lb:\fIt\fR = 0.5

0.45 0.75
0.40 0.69

0.45 0.75
0.25 0.12

0.45 0.73 5lb:\fIt\fR = 1
@
LEGEND=$TMP/legend.dat
draw my.pdf PDE_LQ.dat PDE_HQ.dat $TMP/sol

こんな感じです．

従来の CONV.hpp とかをコピーして, 何かこう基底ファンダメンタルな感じの名前にしておきます. で，改造

// pSolver1DBase.hpp    Solve some equation using some method. Periodic problem
#pragma once                 // この辺は同じ・・・てゆうか，基本，同じ
#include <iostream>
#include <filesystem>
#include <cmath>
#include <vector>
#include <fstream>
#include "pvector.hpp"     
template <typename T>
class pSolver1DBase
{
   protected:
   double xmax = 1.0;
   double dt, CFL;
   double time;
   int nt = 0; // 時間ステップのカウンタ
   class Point
   {
   public:
       Point *pre = nullptr; // ポインタを追加して前のポイントへのリンクを作成
       double X;
       pvector<T> U; // U_n
       Point() : U(3) {}  // デフォルトコンストラクタでサイズ3のpvectorを初期化
       auto &operator[](int index) { return U[index]; }
       const auto &operator[](int index) const { return U[index]; }
   };
   pvector<Point> Data;
   std::ofstream ofs;
public:
   pSolver1DBase(const size_t &nx, const double &dt, const std::function<double(double)> &U0) : dt(dt)
   {
       double C = 1.0;        // 速度
       time = 0.0;            // 初期時間
       double dx = xmax / nx; // 空間刻み幅の計算
       CFL = C * dt / dx;     // CFL条件の計算
       Data.resize(nx);       // 空間分割数に基づいて初期化
       for (size_t i = 0; auto &P : Data)
       {
           P.pre = &Data[i - 1];
           P[nt] = P[nt - 1] = P[nt - 2] = U0(P.X = i * dx); // 初期条件の適用
           i++;
       }
       std::cout << "Booting pSolver1DBaase with nx=" << Data.size() << ", dt=" << dt << std::endl;
   }
   std::string what() const
   {
       return "pSolver1DBase: nx=" + std::to_string(Data.size()) + ", dt=" + std::to_string(dt) +
              ", CFL=" + std::to_string(CFL) + ", time=" + std::to_string(time);
   }
   bool Write(const std::string &filename)
   {
       if (ofs)
           ofs.close(); // 既存のファイルストリームを閉じる
       std::filesystem::path datafile = getenv("HOME");
       datafile /= filename;
       std::filesystem::create_directories(datafile.parent_path());
       ofs.open(datafile);
       if (!ofs)
           throw std::ios_base::failure("Failed to open file: " + datafile.string());
       return Write();
   }
   bool Write()
   {
       try
       {
           ofs << std::endl
               << "# nt= " << nt << " time= " << time << " CFL= " << CFL << std::endl;
           for (auto &v : Data)
               ofs << v.X << " " << v[nt] << std::endl;
           ofs << Data.front().X + xmax << " " << Data.front()[nt] << std::endl; // 最初のデータを追加
       }
       catch (const std::exception &e)
       {
           std::cerr << "Error writing to file: " << e.what() << std::endl;
           return false;
       }
       std::cout << "time[" << time << "]written to file successfully." << std::endl;
       return true;
   }
   virtual void Step(double U_W=0.0)=0;
   virtual void Initialize(void *parm=nullptr) {};
   double get_time() const { return time; }
};

ほとんど元のままですが，

• クラス名は変えたよ. 起動の様子とかみたかったんで，こことかここは，名前を入れたらAIがこうしろっていう
• ここはテンプレートを学んだので，いたずら心でやってみた．
• ここはなんだ？以前は何も書いてなかったので黙ってprivate変数でした．で，private変数は本当にプライベートで，派生型からもアクセスできません．派生型からは見えるけど，外部からは見えない変数が protected 変数です．
• virtual というのをつけると, 「派生型かどうか調べてから，合致するのを使う」となります． virtualでない関数は，定義されていればそのインスタンスの宣言文に多応するクラスの関数が呼び出されます．
• で, virtrual 関数にデフォルトの関数を定義しても良い. ここでは，関数定義{.....} を，ガバッと =0 と書いています．これは
  • 基底クラスをインスタンスにできない
  • 派生クラスを作成する場合, =0 にしてある関数を，必ず絶対抜かりなく定義しなければならぬ．=0が残っている場合，インスタンスにできない.
  • これで，定義し忘れミスを防ぐことができる．
• まあ今のところ使い道がありませんが，こういう，追加の作業が記述できる関数を用意しておくのもありですね．

では，派生してみましょう．簡単のため，最初に作った低精度スキームを書いてみます．

// pSolverLD.hpp
#include "pSolver1DBase.hpp"
class pSolverLD : public pSolver1DBase<double>
{
public:
   pSolverLD(const size_t &nx, const double &dt, const std::function<double(double)> &U0 = [](double x){ return 0.; })
　　　　 : pSolver1DBase(nx, dt, U0)
   {
       std::cout << "Booting pSolverLD..." << std::endl;
   }
   void Step(double U_W = 0.0) override
   {
       for (auto &P : Data) P[nt+1] = P[nt] - CFL * (P[nt] - (*P.pre)[nt]);
       time = (++nt) * dt;
   }
};

まずここで規定クラスを読み込みます．で，派生クラスの定義では，この辺に public と書いてから，基底クラス名を入れます(今回はたまたま基底クラスがテンプレートなのでdoubleで指定)．で，コンストラクトする際に，基底クラスのコンストラクタを呼び出しましょう！で，基底クラスで＝0になっていた関数を，全て定義していきます！このoverride ってのは, 最近は書いた方が絶対良い（省いてもエラーではない）らしいです．

これで，「基本的には基底クラスなんじゃけど」「基底クラスでやること全部やってから」「特殊なのはここね♪」と，書けました.

main()では，まあ普通な感じで使えばよろしい . 読み込むのは派生クラスの定義だけでいいよ:

// main.cpp
#include <iostream>
#include "pSolverLD.hpp"
int main()
{
   size_t nx=10; // 空間分割数
   double dt=0.02; // 時間刻み幅
   size_t ntmax=50; // 最大時間ステップ数
   size_t ntsave=25; // データ保存間隔
   pSolverLD mySolver(nx, dt, [](double x)
                      { return std::pow(std::sin(2.0 * M_PI * x), 5); }); // 初期条件として正弦波を設定
   mySolver.Write("Data/test03/PSOLVER_LD.dat"); // 初期状態の保存
   for(size_t nt=0;nt<ntmax;nt++)
   {
       mySolver.Step(); // 時間ステップを実行
       if ((nt+1) % ntsave == 0) mySolver.Write();
   }
   return 0;
}

実行してみると, 

Booting pSolver1DBaase with nx=10, dt=0.02  ←基底コンストラクタが走って
Booting pSolverLD...　　　　　　　　　　　　←派生コンストラクタが走って
time[0]written to file successfully.        ←Write()は規定のやつ.
time[0.5]written to file successfully.        メッセージは出てないけど, Step()は派生のやつが動く
time[1]written to file successfully.

これで，別のスキーム試したい時には，派生クラスを一つ作って，違う部分だけを書く．ができますね．

draw.sh のスクリプト出しときましょう

#!/bin/bash
#   PDEh.sh: 偏微分方程式I: 周期解
. draw.sh
. standards.sh
COLOR=1
X_AXIS='\fIX\sb\fR1\eb'
X_RANGE='0 1 0.25'
Y_AXIS='\fIu'
Y_RANGE='-1.25 1.25 0.5'
LINE=( Dash:Red:0.003 Solid:Blue:0.002 Solid:Green:0.003 Solid:Cyan:0.003 )
SYMBOL=( Omit )
OPT="--pen-colors=1=DarkRed:2=MediumSeaGreen:3=Blue:4=Magenta:5=GoldenRod3 --symbol-line-width 0.0020"
    ←これが, 色味を変えるコマンド. 色の名前はX11Colorですよ
#解析解のデータが，欲しい
seq 0 0|\
   gawk '
BEGIN{pi=4*atan2(1,1);}
{
time=$1*0.5;
printf("#time= %f\n",time);
for(i=0;i<101;i++)
{
   x=i/100;y=sin(2*pi*(x+time))^5;printf("%f %f\n",x,y);
}
printf("\n");
}
END{printf("\n");}' > $TMP/sol
cat > $TMP/legend.dat <<@
#m=Solid:Black,S=Omit
0.45 1. l:\fIt\fR = 1

0.65 1. l:CFL = 0.2
0.45 0.8    l:\fIi\fR-Sakurai scheme
#m=Dash:Red,S=Omit
0.1 -0.5
0.2 -0.5    l: theory
#m=Solid:Blue,S=Omit
0.1 -0.7
0.2 -0.7    l: \fIN\sbX\eb\fR\r8=\r820
#m=Solid:Green,S=Omit
0.1 -0.9
0.2 -0.9    l: 40
#m=Solid:Cyan,S=Omit
0.1 -1.1
0.2 -1.1    l: 80
@
LEGEND=$TMP/legend.dat
draw HD.pdf $TMP/sol PSOLVER_HD20.dat PSOLVER_HD80.dat PSOLVER_HD320.dat

こんな感じかな.

物理では, \(u\)の総量は，系の総質量とか全エネルギーに対応することが多いので, \(u\)の総量が保存されることは重要です．方程式を見ても(境界周りは別の処理が必要ですが例えば周期領域\(0<x<L\)内の問題だと)\[\int_{0<x<L}\left(\dfrac{\partial u}{\partial t}+c\dfrac{\partial u}{\partial x}\right)\mathrm{d}x=\left(\dfrac{\mathrm{d}}{\mathrm{d}t}\int_{0<x<L}u\mathrm{d}x\right)+c\left[u\right]_{x=0}^{x=L}=0\]で変化しませんが，これは，有限体積法では, \[\int u\left(\Delta_{t},x\right)\mathrm{d} x=\Delta_{x}\sum_{i}\bar{u}_{i}\left(\Delta_{t}\right)=\Delta_{x}\sum_{i}\bar{u}_{i}\left(0\right)+C\sum_{i}\left(F_{i+1/2}-F_{i-1/2}\right)=\Delta_{x}\sum_{i}\bar{u}_{i}\left(0\right)=\int u\left(0,x\right)\mathrm{d} x\]のように\(i\)番目の右数値フラックスと\(i+1\)番目の左数値フラックスがキャンセルして, \(u\)の総量が原理的に保存する，という利点があります．

格子とセルは1個数が違うので, オリジナルの基底クラスの，そこを訂正しましょう.

// pSolverFVBase.hpp    Periodic soltuion by FV method using Cellwise Linear Reconsruction
#pragma once
#include "pSolver1DBase.hpp"
class FVKit
{
public:
   double U;       //セル内平均
   double dU;      //セル内勾配 ∂U/∂x
   double RFlux;   //右境界フラックス/c
   FVKit &operator=(const double &x)
   {
       U = x;
       RFlux = dU = 0;
       return *this;
   }
   operator double() const { return U; }
};
class pSolverFVBase : public pSolver1DBase<FVKit>
{
protected:
   void Grad2Flux()    // Calculate Numerical Flux by Linear Reconstruction
   {
       for (auto &P : Data) P[nt].RFlux=P[nt].U+0.5*(1.-CFL)*P[nt].dU*dx;
   }
public:
   pSolverFVBase(const size_t &nx, const double &dt, const std::function<double(double)> &U0) : pSolver1DBase(nx, dt, U0)
   {
       std::cout << "Booting pSolverFVBase..." << std::endl;
       //  格子点ではなくセルなので1個少ない
       Data.resize(Data.size()-1);
       Data.front().pre=&Data.back();
       Data.back().next=&Data.front();
       for (auto &P : Data) P.X += 0.5 * dx;   //  Xをせる重心に変更
   }
   void Initialize(void *parm=nullptr) override
   {
       std::cout << "pSolverFVBase::Initialize()" << std::endl;
   }
   virtual void getFlux()=0;
   void Step(double U_W = 0.0) override
   {
       getFlux();
       for (auto &P : Data)
           P[nt+1].U = P[nt].U - CFL * (P[nt].RFlux - (*P.pre)[nt].RFlux);
       time = (++nt) * dt;
   }
};

差分だとUだけですが, フラックスとか, 傾きを保存できるようにパワーアップしましょう. で, ここは共通の(4)式を書いておきました．で, データを一個減らして，左右のアクセスを訂正しておきます．なお, back()は配列の最後の要素の参照，front()は配列の最後の要素の参照なので，こうやって繋いでおくと周期配列になります．

数値フラックスの計算は，傾きの計算に依存するので，このクラスを派生して実装しましょう．このへんは, (3)式を実装しています．で，これを実行する前に，ここでgetFlux()関数で傾きを計算しておけばOK,という寸法だ．

例として, MinModのやつを示しておきましょう:

// pSolverFVMinMod.hpp    TVD periodic Finite Volume with MinMod Gradient
#pragma once
#include "pSolverFVBase.hpp"
class pSolverFVMinMod : public pSolverFVBase
{
   double MinMod(const double a,const double b)
   {
       if (a*b<0) return 0.;
       if (std::abs(a)<std::abs(b)) return a;
       return b;
   }
public:
   pSolverFVMinMod(const size_t &nx, const double &dt, const std::function<double(double)> &U0 = [](double x)
                                                   { return 0.; }) : pSolverFVBase(nx, dt, U0)
   {
       std::cout << "Booting pSolverFVMinMod..." << std::endl;
   }
   virtual void getFlux() override
   {
       for (auto &P : Data)
       {
           auto& L=*P.pre;
           auto& R=*P.next;
           P[nt].dU = MinMod(P[nt].U - L[nt].U,R[nt].U-P[nt].U)/dx;
       }
       Grad2Flux();
   }
};

もはや説明するC++文法はないな. まあここくらいか．これは派生クラスで定義されてないので，pSolverFVBaseクラスで定義したやつが実行される．

まあmain()の例

// main.cpp
#include <iostream>
#include "pSolverFVMinMod.hpp"  //MinMod-Limitter
int main()
{
   size_t nx = 20;          // 空間分割数 100(text), 40, 20, 8
   double dt = 0.8 / nx;    // 時間刻み幅 CFL=0.8 に固定しとこ. つまり 0.2/40
   size_t ntmax = 5. / dt;  // 最大時間ステップ数
   size_t ntsave = 5. / dt; // データ保存間隔
   auto I_SIN = [](double x)
   { return std::pow(std::sin(2.0 * M_PI * x), 5); }; // 初期条件として正弦波^5を設定
   auto I_LEVEQUE = [](double x)
   {
       double b=200;
       double xc=x-std::floor(x+0.2)-0.3;
       double y=exp(-b*xc*xc);
       if ((xc>0.3)&&(xc<0.5)) y+=1;
       return y;
   }; // Lebequeのテキストの初期関数
   auto I_DC = [](double x)
   { return (x > 0.8) ? 0 : std::max(2. * x - 0.6, 0.); };
   //--------------------------------------------------
   std::string filename = "Data/test03/TEXT_FVMM20.dat";
   pSolverFVMinMod mySolver(nx, dt, I_LEVEQUE);
   //--------------------------------------------------
   mySolver.Initialize();
   for (size_t nt = 0; nt < ntmax; nt++)
   {
       mySolver.Step(); // 時間ステップを実行
       if ((nt + 1) % ntsave == 0)
           mySolver.Write(filename);
   }
   return 0;
}

ここで定義しているラムダ式は, Levequeの教科書(うちの大学なら無料！)に記載されている例である.