#include "glibw32.h"
#include<stdio.h>
#include<math.h>
#define haikei BLACK
#define L 1000 //x方向の総分割数
const int TOTALTIME=1500; //総時間ステップ数
const double MAXamplitude=1; //波の最大振幅
main(void)
{
int M=TOTALTIME;
double u[L+2],preu[L+2],nextu[L+2],
dt=0.01, dx=0.01, v=1.0, max=MAXamplitude;
ginit(600,200,haikei);
GRAPH gra; gra.view(0,0,600,200);
gra.window(-dx*L/2,-max,dx*L/2,max);
///////初期条件の設定////////////////////////////
//
for(int l=0;l<=L+1;l++)
{
u[l]=exp(-pow((l-L/2)*dx,2));
preu[l]=u[l-v*dt/dx];
}
/////////数値シュミレーション/////////////////////
//
for(int i=1;i<=M;i++)
{
for(int l=1;l<=L;l++)//時間発展
nextu[l]=pow(v*dt/dx,2)*(u[l+1]-2*u[l]+u[l-1])+2*u[l]-preu[l];
nextu[0]=nextu[L]; nextu[L+1]=nextu[1]; //周期境界条件
//
/////////描画部//////////////////////////////
if(i%3==0)
{
gra.frectangle(-dx*L/2,-max,dx*L/2,max,haikei);//消去
gra.moveto(0,-max); gra.lineto(0,max); //縦軸
gra.moveto(-dx*L/2,0);gra.lineto(dx*L/2,0); //横軸
}
gra.moveto(-dx*L/2,u[0]);
for(int l=1;l<=L+1;l+=3)
{
gra.lineto(dx*(l-L/2),u[l],GREEN);
}
gra.moveto(-dx*L/2,u[0]);

/////////////入れ替え preu <- u, u <- nextu ///////////////
for(int l=0;l<L+2;l++)
{
preu[l]=u[l];u[l]=nextu[l];
}
}
gend();
}