数値計算法レポ2ルンゲクッタ版(OK)

「数値計算法レポ2ルンゲクッタ版(OK)」の編集履歴(バックアップ)一覧はこちら

数値計算法レポ2ルンゲクッタ版(OK)」(2011/07/08 (金) 10:57:07) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

-正直,[[NSMRのやつ>http://www47.atwiki.jp/cscd/pages/43.html]]のほうがいいプログラムだと思う。配列使ってないという意味で。 -[[ファイル名の自動化]]を使用している。 -後半はgnuplot用のマクロファイルを作っていて,それをgvで表示するところまで入れてある。 #highlight(linenumber,c){{ #include<stdio.h> #include<stdlib.h> #include<math.h> int i,j; double x[9999],y[9999]; double kx11[9999],kx12[9999],kx13[9999],kx14[9999]; double ky11[9999],ky12[9999],ky13[9999],ky14[9999]; double kx21[9999],kx22[9999],kx23[9999],kx24[9999]; double ky21[9999],ky22[9999],ky23[9999],ky24[9999]; double vx[9999],vy[9999]; double r,dt,p; FILE *fpOUT; FILE *fpPLT; double Z;//ターゲットの原子番号 double m,a,c,hbarc;//入射粒子の質量,微細構造定数,hバーかける光速 double A; double ax,ay; // http://questionbox.jp.msn.com/qa1138947.html FILE *my_fopenw(int x) { // int x; if (x >= 0 /* && x < 2 */) { char file[256]; sprintf(file, "./RK%d.dat",x); return fopen(file, "w"); } return (FILE *)NULL; } FILE *my_fopenr(int x) { // int x; if (x >= 0 /* && x < 2 */) { char file[256]; sprintf(file, "./RK%d.dat",x); return fopen(file, "r"); } return (FILE *)NULL; } int main(){ //パラメーター,定数 はじまり dt = 1.0;//zs c = 3.0*100.0; // speed of light 3E08 fm/zs p = 5.3;//MeV/c a = 1.0/137.0;//fine constructure constant hbarc = 197.0; //[MeV*fm]=197MeV*fm m = 3.727*1000.0; //α粒子の質量=3.727379109(93) GeV/c2 Z = 79.0; //for Au nuclear A = Z*a*hbarc/m; x[0] = -2000.0;//fm vx[0]=p*c/m; vy[0]=0; j=0; //おわり // 衝突係数の入力 printf("Enter the impact parameter (fm):\n"); scanf("%lf",&y[0]); // ファイル名の設定↓ while((fpOUT = my_fopenr(j)) != NULL){ j++; } fpOUT = my_fopenw(j); printf("Create new file RK%d.dat\n",j); // ファイル名の設定↑ fprintf(fpOUT,"%lf %lf\n",x[0],y[0]); for(i = 1;i <= 9999;i++){ r = sqrt(x[i-1]*x[i-1]+y[i-1]*y[i-1]); kx11[i-1] = A*x[i-1]/(r*r*r); kx21[i-1] = vx[i-1]; kx12[i-1] = A*(x[i-1]+kx21[i-1]*dt/2)/(r*r*r); kx22[i-1] = vx[i-1]+kx11[i-1]*dt/2; kx13[i-1] = A*(x[i-1]+kx22[i-1]*dt/2)/(r*r*r); kx23[i-1] = vx[i-1]+kx12[i-1]*dt/2; kx14[i-1] = A*(x[i-1]+kx23[i-1]*dt)/(r*r*r); kx24[i-1] = vx[i-1]+kx12[i-1]*dt; vx[i] = (kx11[i-1]+2*kx12[i-1]+2*kx13[i-1]+kx14[i-1])*dt/6 + vx[i-1]; x[i] = (kx21[i-1]+2*kx22[i-1]+2*kx23[i-1]+kx24[i-1])*dt/6 + x[i-1]; ky11[i-1] = A*y[i-1]/(r*r*r); ky21[i-1] = vy[i-1]; ky12[i-1] = A*(y[i-1]+ky21[i-1]*dt/2)/(r*r*r); ky22[i-1] = vy[i-1]+ky11[i-1]*dt/2; ky13[i-1] = A*(y[i-1]+ky22[i-1]*dt/2)/(r*r*r); ky23[i-1] = vy[i-1]+ky12[i-1]*dt/2; ky14[i-1] = A*(y[i-1]+ky23[i-1]*dt)/(r*r*r); ky24[i-1] = vy[i-1]+ky12[i-1]*dt; vy[i] = (ky11[i-1]+2*ky12[i-1]+2*ky13[i-1]+ky14[i-1])*dt/6 + vy[i-1]; y[i] = (ky21[i-1]+2*ky22[i-1]+2*ky23[i-1]+ky24[i-1])*dt/6 + y[i-1]; fprintf(fpOUT,"%lf %lf\n",x[i],y[i]); } fclose(fpOUT); // Plot with GNUPLOT fpPLT = fopen("RK.plt", "w"); fprintf(fpPLT,"set term post eps colour\n"); fprintf(fpPLT,"set output \"RK.eps\"\n"); fprintf(fpPLT,"set xlabel \"x [fm]\"\n"); fprintf(fpPLT,"set ylabel \"y [fm]\"\n"); fprintf(fpPLT,"set key below\n"); fprintf(fpPLT,"plot [-2000:2000] [-2000:2000]\"RK0.dat\" w l"); for(i=1;i<=j;i++){ fprintf(fpPLT,", \"RK%d.dat\" w l", i); } fprintf(fpPLT,"\n"); fprintf(fpPLT,"set term wxt"); fclose(fpPLT); system("gnuplot RK.plt"); system("gv RK.eps"); } }} #image(RK.gif,title=result) コメント募集中 #comment_num2()
-正直,[[NSMRのやつ>http://www47.atwiki.jp/cscd/pages/43.html]]のほうがいいプログラムだと思う。配列使ってないという意味で。 -[[ファイル名の自動化]]を使用している。 -後半はgnuplot用のマクロファイルを作っていて,それをgvで表示するところまで入れてある。 #highlight(linenumber,c){{ #include<stdio.h> #include<stdlib.h> #include<math.h> int i,j; double x[9999],y[9999]; double kx11[9999],kx12[9999],kx13[9999],kx14[9999]; double ky11[9999],ky12[9999],ky13[9999],ky14[9999]; double kx21[9999],kx22[9999],kx23[9999],kx24[9999]; double ky21[9999],ky22[9999],ky23[9999],ky24[9999]; double vx[9999],vy[9999]; double r,dt,p; FILE *fpOUT; FILE *fpPLT; double Z;//ターゲットの原子番号 double m,a,c,hbarc;//入射粒子の質量,微細構造定数,hバーかける光速 double A; double ax,ay; // http://questionbox.jp.msn.com/qa1138947.html FILE *my_fopenw(int x) { // int x; if (x >= 0 /* && x < 2 */) { char file[256]; sprintf(file, "./RK%d.dat",x); return fopen(file, "w"); } return (FILE *)NULL; } FILE *my_fopenr(int x) { // int x; if (x >= 0 /* && x < 2 */) { char file[256]; sprintf(file, "./RK%d.dat",x); return fopen(file, "r"); } return (FILE *)NULL; } int main(){ //パラメーター,定数 はじまり dt = 1.0;//zs c = 3.0*100.0; // speed of light 3E08 fm/zs p = 5.3;//MeV/c a = 1.0/137.0;//fine constructure constant hbarc = 197.0; //[MeV*fm]=197MeV*fm m = 3.727*1000.0; //α粒子の質量=3.727379109(93) GeV/c2 Z = 79.0; //for Au nuclear A = Z*a*hbarc/m; x[0] = -2000.0;//fm vx[0]=p*c/m; vy[0]=0; j=0; //おわり // 衝突係数の入力 printf("Enter the impact parameter (fm):\n"); scanf("%lf",&y[0]); // ファイル名の設定↓ while((fpOUT = my_fopenr(j)) != NULL){ j++; } fpOUT = my_fopenw(j); printf("Create new file RK%d.dat\n",j); // ファイル名の設定↑ fprintf(fpOUT,"%lf %lf\n",x[0],y[0]); for(i = 1;i <= 9999;i++){ r = sqrt(x[i-1]*x[i-1]+y[i-1]*y[i-1]); kx11[i-1] = A*x[i-1]/(r*r*r); kx21[i-1] = vx[i-1]; kx12[i-1] = A*(x[i-1]+kx21[i-1]*dt/2)/(r*r*r); kx22[i-1] = vx[i-1]+kx11[i-1]*dt/2; kx13[i-1] = A*(x[i-1]+kx22[i-1]*dt/2)/(r*r*r); kx23[i-1] = vx[i-1]+kx12[i-1]*dt/2; kx14[i-1] = A*(x[i-1]+kx23[i-1]*dt)/(r*r*r); kx24[i-1] = vx[i-1]+kx12[i-1]*dt; vx[i] = (kx11[i-1]+2*kx12[i-1]+2*kx13[i-1]+kx14[i-1])*dt/6 + vx[i-1]; x[i] = (kx21[i-1]+2*kx22[i-1]+2*kx23[i-1]+kx24[i-1])*dt/6 + x[i-1]; ky11[i-1] = A*y[i-1]/(r*r*r); ky21[i-1] = vy[i-1]; ky12[i-1] = A*(y[i-1]+ky21[i-1]*dt/2)/(r*r*r); ky22[i-1] = vy[i-1]+ky11[i-1]*dt/2; ky13[i-1] = A*(y[i-1]+ky22[i-1]*dt/2)/(r*r*r); ky23[i-1] = vy[i-1]+ky12[i-1]*dt/2; ky14[i-1] = A*(y[i-1]+ky23[i-1]*dt)/(r*r*r); ky24[i-1] = vy[i-1]+ky12[i-1]*dt; vy[i] = (ky11[i-1]+2*ky12[i-1]+2*ky13[i-1]+ky14[i-1])*dt/6 + vy[i-1]; y[i] = (ky21[i-1]+2*ky22[i-1]+2*ky23[i-1]+ky24[i-1])*dt/6 + y[i-1]; fprintf(fpOUT,"%lf %lf\n",x[i],y[i]); } fclose(fpOUT); // Plot with GNUPLOT fpPLT = fopen("RK.plt", "w"); fprintf(fpPLT,"set term post eps colour\n"); fprintf(fpPLT,"set output \"RK.eps\"\n"); fprintf(fpPLT,"set xlabel \"x [fm]\"\n"); fprintf(fpPLT,"set ylabel \"y [fm]\"\n"); fprintf(fpPLT,"set key below\n"); fprintf(fpPLT,"plot [-2000:2000] [-2000:2000]\"RK0.dat\" w l"); for(i=1;i<=j;i++){ fprintf(fpPLT,", \"RK%d.dat\" w l", i); } fprintf(fpPLT,"\n"); fprintf(fpPLT,"set term wxt"); fclose(fpPLT); system("gnuplot RK.plt"); system("gv RK.eps"); } }} 実行結果の例 #image(RK.gif,title=result) コメント募集中 #comment_num2()

表示オプション

横に並べて表示:
変化行の前後のみ表示: