To simulate the changes of the distribution of the substance concentration with the passage of time , we recurrently calculate the status of only time of ⊿t after from the current state. This method is called difference approximation. We will solve the equation obtained by it. "Explicit" means the method of determining the next state from the current state. There is also implicit method and it has advantage to be able to take large time increment ⊿t. However, it is complex to solve equations for each time. Compared to implicit method you will be able to simulate the diffusion in a simple program if you use the explicit method.
This source code is protected by the GNU GPL license. Please visit here for the terms of the GNU GPL license.
//
// 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;
}