プログラム例
#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
こんなかんじで.