「プログラムソース」
９−１、-------------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#define D 0.01

unsigned rndnum=13;
unsigned irnd(void);

int main(){
	
	int counter,trial_times,randum_number;
	double x,y,z=0;
	double input_x,input_y,into_z;

	printf("値を入力して下さい。(0<=x,y<=1)\n");
	printf("x=");
	scanf("%lf",&x);
	printf("y=");
	scanf("%lf",&y);
	printf("試行回数を入力して下さい。\n");
	printf("trialtimes=");
	scanf("%d",&trial_times);
	printf("計算中・・・");

	input_x=x;
	input_y=y;
	
	for(counter=1;counter<=trial_times;counter++){
		x=input_x;
		y=input_y;
		while(((0<=x && x<=1) && 0<=y) && y<=1){
			randum_number=(int)(14*(irnd()/32767.1));
			if(randum_number==0||randum_number==1||randum_number==2||randum_number==3||randum_number==4)
				x+=D;
			else if(randum_number==5||randum_number==6||randum_number==7||randum_number==8||randum_number==9)
				x-=D;
			else if(randum_number==10||randum_number==11)
				y+=D;
			else if(randum_number==12||randum_number==13)
				y-=D;
			else{
				printf("エラー1\n");
				break;
			}
		}
		
		if(x>1)
			z=y;
		else if(x<0)
			z=1-y;
		else if(y>1)
			z=x;
		else if(y<0)
			z=1-x;
		else{
			printf("エラー2\n");
			break;
		}
		into_z+=z;
		
	}
	
	printf("f(%f,%f)=%f\n",input_x,input_y,into_z/trial_times);
	return 0;
}

unsigned irnd()
{
	rndnum=(rndnum*109+1021)%32768;
	return(rndnum);
}

--------------------------------------------------------------------------------

「実行結果」
９−１--------------------------------------------------------------------------


C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.5
y=0.5
試行回数を入力して下さい。
trialtimes=10
計算中・・・f(0.500000,0.500000)=0.369000

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.5
y=0.5
試行回数を入力して下さい。
trialtimes=100
計算中・・・f(0.500000,0.500000)=0.402400

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.5
y=0.5
試行回数を入力して下さい。
trialtimes=10000
計算中・・・f(0.500000,0.500000)=0.404974

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.3
y=0.5
試行回数を入力して下さい。
trialtimes=10
計算中・・・f(0.300000,0.500000)=0.432000

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.3
y=0.5
試行回数を入力して下さい。
trialtimes=100
計算中・・・f(0.300000,0.500000)=0.422300

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.3
y=0.5
試行回数を入力して下さい。
trialtimes=10000
計算中・・・f(0.300000,0.500000)=0.420037

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.5
y=0.3
試行回数を入力して下さい。
trialtimes=10
計算中・・・f(0.500000,0.300000)=0.465000

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.5
y=0.3
試行回数を入力して下さい。
trialtimes=100
計算中・・・f(0.500000,0.300000)=0.449700

C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\試験ソース>monntekaruro2
値を入力して下さい。(0<=x,y<=1)
x=0.5
y=0.3
試行回数を入力して下さい。
trialtimes=10000
計算中・・・f(0.500000,0.300000)=0.448017


--------------------------------------------------------------------------------

「考察」
与えられた偏微分方程式をモンテカルロ法を用いて解くプログラムを作成した。

入力されたx,yの値におけるf(x,y)を求めるには、そのx,yを格子点における始点と考えて、出力される乱数にしたがって移動する試行を繰り返す。
境界に達する時の値を複数回の試行によって導きだし、それらの平均値をとる事で収束する値を求める。

ここでは、偏微分方程式の係数に注目して両辺を14で割り、14で乱数を割った時のあまりを利用して格子点の移動を行う。
ただし、係数からこの場合はxは5/14の確率で一方向に移動し、yは1/7の確率で移動する。

ソース内では発生させる乱数をrandum_numerとして、その値にしたがって格子点の移動を行っている。

境界条件に格子点が達した場合は、その値をinto_zに次々と加えていき、最後に平均値をとるためにinto_zに入れたループ回数で割っている。


--------------------------------------------------------------------------------