メインコンテンツに移動

解答例1

プログラム例

#include <array>
#include <vector>
#include <boost/filesystem.hpp>
class txData
{
private:
  static const int step=2;
  std::array<std::vector<double>,step> data;
  int nx;
public:
  std::vector<double>& operator[](int nt)
  {
    return data[nt%step];
  }
  txData(int N):nx(N)
  {
    for(int nt=0;nt<step;nt++) data[nt].resize(N);
  }
  int size()
  {
    return nx;
  }
};
class Solver
{
  txData U,V;
  double c=1.0,cfl;
  double dt,dx;
  std::vector<double> X;
  std::ofstream ofs;
public:
  Solver(double xmax,int nx,const double dt,std::string init,std::string filename)
  :U(nx+1),V(nx+1),dt(dt)
  {
    if (nx%2!=0) throw(std::runtime_error("nxが偶数じゃないぞタコ"));
    boost::filesystem::path folder=getenv("HOME");
    boost::filesystem::current_path(folder/"tes_test");
    dx=xmax/nx; //0,1,2,...,nx
    X.resize(U.size());
    for(int i=0;i<X.size();i++) X[i]=-xmax/2+dx*i;
    cfl=c*dt/dx;
    if (cfl>1.0) throw(std::runtime_error("cflでかすぎぢゃボケ"));
    //色んな初期値で遊ぶ
    if (init=="kaku")
    {//カクカクした初期値の例
      for(int i=0;i<X.size();i++)
      {
        U[0][i]=(std::abs(X[i])<=1.0)?1.0:0.0;
        V[0][i]=(std::abs(X[i])<=1.0)?0.0:0.0;
      }
    }
    else
    {//ぬるぬるした初期値の例
      for(int i=0;i<X.size();i++)
      {
        U[0][i]=std::exp(-5*X[i]*X[i]);
        V[0][i]=0.0;
      }
    }
    ofs.open(filename);
  }
  void Step(int nt)
  {
    //境界条件を設定
     U[nt+1][0]=0.;
     V[nt+1][X.size()]=0.;
    //解きさらせ
    for(int i=1;i<=X.size();i++)
      U[nt+1][i]=U[nt][i]-cfl*(V[nt][i]-V[nt][i-1]);
    for(int i=0;i<X.size();i++)
      V[nt+1][i]=V[nt][i]-cfl*(U[nt][i+1]-U[nt][i]);
  }
  void Write(int nt)
  {
    for(int i=0;i<X.size();i++)
      ofs << X[i] << " " << U[nt][i] << std::endl;
    ofs<<std::endl;
  }
};

ほんでmain()では

double xmax=....
double dt=...
Solver mySolver(xmax,nx,dt,"kaku","my_wave.data");
mySolver.Write(0);
for(int nt=0;nt<ntmax;nt++)
{
  mySolver.Step(nt);
  if ((nt+1)%ntsave==0) mySolver.Write(nt+1);
}

こんな感じですね.グラフ描きは

graph -Tps -y -0.5 1.5 -x -4 +4 -C my_wave.data |ps2pdf -dEPSCrop - new.pdf

こんなかんじで.