2010夏合宿山田 - P.S.P.T. 計算物理研究会

ローレンツカオス

• 発表者 28期　山田
• 2010年度夏合宿
• 使用言語 C言語
• 添付資料

使用したパワーポイント

• 添付ソースコード

ソースファイル_LorenzAttractor_c.zip

概要 結果 考察 感想 参考文献 使用方法 ソースコード

---

概要

ローレンツ方程式とは

x : 対流運動の強さに比例した変数
y : 上昇流と下降流の温度差に比例した変数
z : 鉛直方向の温度変化が1 次関数からどの程度ずれているかを表す変数

プラントル数
規格化されたレーリイ数
与えられた容器に関係する定数

結果

下の図は
= 10.0
= 28.0
= 8.0/3.0
であり、初期値が
赤線が
(0)=0.0
(0)=1.0
(0)=1.0
緑線が
(0)=0.0
(0)=1.000001
(0)=1.0
の相空間をプロットしたものです。

そして、次の図はとをそれぞれtについてプロットしたものです。

t=30近辺から大きくズレていってるのがわかります。

考察

Gnuplotで描いた図からはローレンツカオスを観測できたと思わえる。

感想

これからは、紙とペンで解析出来るように勉強していきたい。

参考文献

---

使用方法

下のサイトを参考にしてGnuplotを直接呼び出せるようにしてください。
物理のかぎしっぽ_グラフ解析・ツール

ソースコード

#include <stdio.h>
 
#define dt 0.005
#define N 2 //初期値のかず
 
double f1(double y1, double y2, double y3, double s);
double f2(double y1, double y2, double y3, double r); 
double f3(double y1, double y2, double y3, double b);
 
int main()
{
//定数
	//①
	double s = 10.0; 
	double r = 28.0;
	double b = 8.0/3.0;
	//②
//	double s = 10.4; 
//	double r = 81.6;
//	double b = 6.8;	
	//③
//	double s = 7.8; 
//	double r = 49.7;
//	double b = 4.4;	
 
	int i;
	int loop_max = 70.0;
	double t;
	double y1[N] , y2[N] , y3[N];
 
        double k1[4], k2[4], k3[4];
 
	FILE *fp;
 
	//初期条件
	t = 0.0;
 
	y1[0] = 0.0;
	y2[0] = 1.0;
	y3[0] = 1.0;
 
	y1[1] = 0.0;
	y2[1] = 1.000001;
	y3[1] = 1.0;
  /*
	y1[2] = 1.0;
	y2[2] = 1.0;
	y3[2] = 1.0;
*/
 	fp = fopen("LorenzAttractor.dat", "w");
 
	for(t=0; t<loop_max; t += dt){
		fprintf(fp,"%f ",t);
		for (i = 0; i < N; i++) {
 
			k1[0] = f1(y1[i], y2[i], y3[i], s)*dt;
			k2[0] = f2(y1[i], y2[i], y3[i], r)*dt;
			k3[0] = f3(y1[i], y2[i], y3[i], b)*dt;
 
			k1[1] = f1(y1[i]+k1[0]/2.0, y2[i]+k2[0]/2.0, y3[i]+k3[0]/2.0, s)*dt;
		        k2[1] = f2(y1[i]+k1[0]/2.0, y2[i]+k2[0]/2.0, y3[i]+k3[0]/2.0, r)*dt;
			k3[1] = f3(y1[i]+k1[0]/2.0, y2[i]+k2[0]/2.0, y3[i]+k3[0]/2.0, b)*dt;
 
			k1[2] = f1(y1[i]+k1[1]/2.0, y2[i]+k2[1]/2.0, y3[i]+k3[1]/2.0, s)*dt;
		        k2[2] = f2(y1[i]+k1[1]/2.0, y2[i]+k2[1]/2.0, y3[i]+k3[1]/2.0, r)*dt;
			k3[2] = f3(y1[i]+k1[1]/2.0, y2[i]+k2[1]/2.0, y3[i]+k3[1]/2.0, b)*dt;
 
			k1[3] = f1(y1[i]+k1[2], y2[i]+k2[2], y3[i]+k3[2], s)*dt;	
			k2[3] = f2(y1[i]+k1[2], y2[i]+k2[2], y3[i]+k3[2], r)*dt;	
			k3[3] = f3(y1[i]+k1[2], y2[i]+k2[2], y3[i]+k3[2], b)*dt;	
 
			y1[i] += (k1[0]+2.0*k1[1]+2.0*k1[2]+k1[3])/6.0; 	
		        y2[i] += (k2[0]+2.0*k2[1]+2.0*k2[2]+k2[3])/6.0;   
			y3[i] += (k3[0]+2.0*k3[1]+2.0*k3[2]+k3[3])/6.0;
 
			fprintf(fp,"%f %f %f ", y1[i], y2[i], y3[i]);	
		}	
		fprintf(fp, "\n");                                                                                        
	}
	fclose(fp);
	system("wgnuplot -persist setting.plt");
	system("wgnuplot -persist setting2.plt");
}                                                
 
double f1(double y1, double y2, double y3, double s)
{
	return s*(-y1+y2);
}                                                                      
double f2(double y1, double y2, double y3, double r)
{
	return r*y1 -y2 -y1*y3;
}
double f3(double y1, double y2, double y3, double b)
{
	return -b*y3 +y1*y2;
}