拡散方程式と陽解法

時間の経過に伴う物質濃度の分布の変化をシミュレーションするために、現在の状態から⊿ｔだけあとの分布の状態を求めます。これには差分近似というものを用い、それによって得られた方程式を解いていきます。現在の状態から次の状態を求める方法が陽解法となります。他に陰解法というものがありますが説明が難しいので省きます。陰解法には時間増分⊿ｔを大きくとれるなどのメリットが有るのですが、各時間ごとに連立方程式を解く必要があり複雑です。それに比べると陽解法を使えば簡単なプログラムで拡散をシミュレートすることができます。

ソースコード

このソースコードはGNU GPLライセンスにより保護されています。GNU GPLライセンスの条項に関してはこちらをご覧ください。

//
//  diffusion
//
//  Copyright 2015 Yuki Mochizuki
//

#include <iostream>
#include <stdio.h>
#include "opencv/cv.h"  
#include "opencv/highgui.h"  

#define NMAX_X 100
#define NMAX_Y 100
#define N 2000000

double r(double,double,double,double);
double f(double,double);
double f2(double,double);
double f3(double,double,int);

int main(){
    cv::VideoWriter writer("kakusan.avi", CV_FOURCC_DEFAULT, 30, cv::Size(NMAX_X, NMAX_Y), true);
    char *output ="noudo.txt";
    FILE* fpo;
    fpo=std::fopen( output, "w" );
    
    cv::Mat img=cv::Mat(cv::Size(NMAX_X,NMAX_Y),CV_8UC3);
    cv::namedWindow("img", CV_WINDOW_AUTOSIZE);
    
    int i,j,k,p;
    double ax,ay,dx,dy,dt,a;
    double U[NMAX_X+1][NMAX_Y+1],UU[NMAX_X+1][NMAX_Y+1];
    double MaxValue;
    
    dt=1.0/1000000.0;
    dx=1.0/(double)NMAX_X;
    dy=1.0/(double)NMAX_Y;
    a=1.0/10.0;
    ax=1.0/dx/dx;
    ay=1.0/dy/dy;
    
    FILE* gp;
    
    if(dt>pow((dx*dy),2)/(dx*dx+dy*dy)/2/a){
        printf("dt is too big to converge\n");
        return -1;
    }
    
    
    for(i=0;i<=NMAX_X;i++){
        for(j=0;j<=NMAX_Y;j++){
            U[i][j]=f3((double)i/(double)NMAX_X,(double)j/(double)NMAX_Y,5);
            UU[i][j]=0.0;
        }
    }
    
    gp = popen("/usr/local/bin/gnuplot -persist","w");
    fprintf(gp, "set terminal x11\n");
    
    
    for(k=0;k<=N;k++){
        
        if(k%500==0){
            img=cv::Scalar(0,0,0);
            MaxValue=0;
            for(i=0;i<=NMAX_X;i++){
                for(j=0;j<=NMAX_Y;j++){
                    if(MaxValue<U[i][j])MaxValue=U[i][j];
                }
            }
            
            for(i=0;i<=NMAX_X;i++){
                for(int j=0;j<NMAX_Y;j++){
                    p = (int)img.step*((int)j)+(((int)i)*3);
                    img.data[p+1] =((int)(U[i][j]*255/MaxValue))%256 ;
                }
            }
            
            cv::imshow("img", img);
            writer << img;
            if(cv::waitKey(1)!=-1){
                fclose(fpo);
                pclose(gp);
                return 0;
            }
            printf("\n");
        }
        
        
  
        for(i=1;i<=NMAX_X-1;i++){
            for(j=1;j<=NMAX_Y-1;j++){
                UU[i][j]=a*dt*(ax*U[i+1][j]+ax*U[i-1][j]+ay*U[i][j+1]+ay*U[i][j-1])+(1.0-2.0*a*dt*(ax+ay))*U[i][j];
            }
        }
        
        for(i=0;i<=NMAX_X;i++){
            for(j=0;j<=NMAX_Y;j++){
                U[i][j]=UU[i][j];
            }
        }
        
        for(i=NMAX_X/2-4;i<=NMAX_X/2+4;i++){
            for(j=NMAX_Y/2-5;j<=NMAX_Y/2+5;j++){
                if(f3((double)i/(double)NMAX_X,(double)j/(double)NMAX_Y,5))U[i][j]=1;
            }
        }
        for(i=0;i<=NMAX_X;i++){
            for(j=0;j<=NMAX_Y;j++){
                if(!f3((double)i/(double)NMAX_X,(double)j/(double)NMAX_Y,50))U[i][j]=0;
            }
        }
        
        
        if(k%10000==0){
            
            for(int ro=0;ro<=NMAX_X;ro+=1){
                    fprintf(fpo, "%d %f\n",ro,U[ro][(int)NMAX_Y/2]);
                printf( "%d %f\n",ro,U[ro][(int)NMAX_Y/2]);
            }
            fprintf(fpo, "\n\n");
            fflush(fpo);
            
            fprintf(gp, "plot ");
            for(int i=0;i<=k/10000;i++){
                fprintf(gp, "\"noudo.txt\" index %d using 1:2 title \"%d\" with lines,",i,i);
                //fprintf(gp, "\"noudo.txt\" index %d using 1:2 title \"%d\" with lp,",i,i);
            }
            fprintf(gp, "\n");
            fflush(gp);
        }
    }
    
    
    fprintf(gp, "\n");
    pclose(gp);
    fclose(fpo);

}

double r(double x1,double y1,double x2,double y2){
    return pow((pow((x1-x2),2)+pow((y1-y2),2)),0.5);
}
double f3(double x,double y,int rr){
    if(r(x*100,y*100,50,50)<rr) return 1;
    else return 0;
}