一次元熱伝導方程式のシミュレーション

heat_eq_1dim.cpp

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/* 1次元熱伝導方程式
 * u_t = ALPHA * u_xx
 * ただし MIN_X <= x <= MAX_X
 * 境界条件
 * u(MIN_X,t) = U_LEFT
 * u(MAX_X,t) = U_RIGHT */

#define N_SPLIT	31	// 分割数-1
#define MIN_X	0.0	// xの最小値
#define MAX_X	1.0	// xの最大値
#define U_LEFT	0.0	// uの左の固定値
#define U_RIGHT	0.0	// uの右の固定値
#define ALPHA	1.0	// 熱伝導率

// t,xの刻み幅
#define dt	1e-4
#define dx ((MAX_X-MIN_X)/N_SPLIT)

// シミュレーションを行う関数(n_stepはすすめるステップ数)
void calc(double *u,int n_step);
// 画面にデータを表示する関数
void disp(double *u,double t);

int main(void){
	// x = MAX_Xを含めるために要素数を一つ増やしている
	double u[N_SPLIT+1];
	
	// 初期条件は自分で決める
	for(int i = 0 ; i <= N_SPLIT ; i++ ){
		double x = (MAX_X - MIN_X)*(double)i / N_SPLIT + MIN_X;
		
		if( x < 0.5 )
			u[i] = 2*x;
		else
			u[i] = 2*(1-x);
	}
	
	for(int i = 0 ; i < 100 ; i++ ){
		disp(u,i*dt);
		calc(u,1);
	}
	
	return 0;
}

void calc(double *u,int n_step){
	// (t+dt)のデータ．一時保存用
	double v[N_SPLIT+1];
	
	// 一応境界条件を入れておく
	u[0] = U_LEFT;
	u[N_SPLIT] = U_RIGHT;
	
	v[0] = U_LEFT;
	v[N_SPLIT] = U_RIGHT;
	
	for(int i = 0 ; i < n_step ; i++ ){
		// 差分法によって得られた式を実行
		// v[0]とv[N_SPLIT]には適用されないので注意
		for(int x = 1 ; x < N_SPLIT ; x++ )
			v[x] = u[x] + ALPHA*(dt/(dx*dx))*(u[x+1]-2*u[x]+u[x-1]);
		
		// uにコピー
		for(int x = 0 ; x <= N_SPLIT ; x++ )
			u[x] = v[x];
	}
}

void disp(double *u,double t){
	printf("%lf,",t);
	for(int i = 0 ; i < N_SPLIT ; i++ )
		printf("%lf,",u[i]);
	printf("%lf\n",u[N_SPLIT]);
}
/*
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

u = pd.read_csv("output.csv",header=None)
x = np.linspace(0,1,num=32)

for i in range(100):
    plt.plot(x,u.ix[i,1:])
    plt.ylim(-1.0,1.0)
    plt.xlim(0.0,1.0)
    
    plt.pause(0.1)
    #plt.savefig(".\\a\\%d.png" % i )
    plt.clf()
*/