pt
// program parabolic1
#include <iostream>
using namespace std;
#include "cpgplot.h"
void initparam(double u[],
int & nx,
double &dx,
double & dt,
int & niter,
int &gmode)
{
cerr << "Enter nx, tmax, dt:";
double tmax;
cin >> nx >> tmax >> dt;
dx = 1.0/nx;
niter = (tmax+dt/2)/dt;
cerr << "Enter graphics mode 0: no G\n"
<< " 1: animation\n"
<< " 2: no-erase)\n";
cin >> gmode;
cout << "nx = " << nx<< " dt = " << dt << " niter = "<< niter<< "\n";
cout << "gmode = " << gmode << "\n";
int i;
for (i=0;i<=nx; i++) u[i] = 0;
u[(nx+1)/2] = nx;
}
void push_system(double u[],
double unew[],
int nx,
double dx,
double dt)
{
double lambda = dt/(dx*dx);
int i;
for (i=1; i<nx; i++){
unew[i] = u[i] + lambda*(u[i-1]-2*u[i]+u[i+1]);
}
for(i=1; i<nx; i++){
u[i] = unew[i];
}
}
void initgraph()
{
if(cpgopen("?") !=1 ) exit(-1);
cpgask(0);
}
void show_graphics(double u[],
int nx,
double dx,
int gmode,
double xmin,
double xmax,
double ymin,
double ymax,
int first)
{
if (gmode == 0) return;
cpgbbuf();
if (first == 1){
cpgpage();
cpgenv(xmin, xmax, ymin, ymax,0,-2);
}
if (gmode == 1) {
cpgeras();
}
cpgbox("bcnst", 0.0, 0, "bcnst", 0.0, 0);
cpglab("X", "Y", " ");
cpgmove(0.0, u[0]);
int i;
for(i=1;i<nx+1;i++){
cpgdraw(i*dx, u[i]);
}
cpgebuf();
}
int main()
{
#define NMAX 10001
static double u[NMAX], unew[NMAX];
double dx, dt;
int nx, niter, gmode, iter;
initparam(u, nx, dx, dt, niter, gmode);
if (gmode != 0) initgraph();
show_graphics(u, nx, dx, gmode, 0.0, 1.0, 0.0, 10.0, 1);
for(iter = 1; iter <=niter; iter++){
push_system(u, unew, nx, dx, dt);
show_graphics(u, nx, dx, gmode,0.0, 1.0, 0.0, 10.0, 0);
}
return 0;
}
これをちょっと動かしてみよう。
コンパイルには
g++ -I/usr/local/include -o cparabolic1 cparabolic1.C \
-L/usr/local/pgplot -lcpgplot -lpgplot \
-L/usr/X11R6/lib -lX11 -lg2c -lm
でいいはずである。 g++ のバージョンによっては、 -lg2c が -lf2c でない
といけないかもしれない。また、 PGPLOT をインストールした場所が違うとそ
のあたりが変わってくる。