version 1.5, 2013/02/23 06:01:15 |
version 1.7, 2013/03/07 03:00:43 |
|
|
/* |
/* |
License: LGPL |
License: LGPL |
Ref: Copied from this11/misc-2011/A1/wishart/Prog |
Ref: Copied from this11/misc-2011/A1/wishart/Prog |
$OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.4 2013/02/20 05:20:49 takayama Exp $ |
$OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.6 2013/03/01 05:26:25 takayama Exp $ |
*/ |
*/ |
#include <stdio.h> |
#include <stdio.h> |
#include <stdlib.h> |
#include <stdlib.h> |
#include <math.h> |
#include <math.h> |
Line 22 static struct SFILE *Gf =NULL; |
|
Line 22 static struct SFILE *Gf =NULL; |
|
static struct SFILE *Df =NULL; |
static struct SFILE *Df =NULL; |
int MH_P95=0; |
int MH_P95=0; |
int MH_Verbose=0; |
int MH_Verbose=0; |
static mypower(int x,int n) { |
static int mypower(int x,int n) { |
int a,i; |
int a,i; |
a = 1; |
a = 1; |
for (i=0; i<n; i++) a = a*x; |
for (i=0; i<n; i++) a = a*x; |
return(a); |
return(a); |
} |
} |
mh_gopen_file() { |
int mh_gopen_file() { |
FILE *fp; |
FILE *fp; |
char fname[1024]; |
char fname[1024]; |
int i; |
int i; |
Line 40 mh_gopen_file() { |
|
Line 40 mh_gopen_file() { |
|
fprintf(stderr,"MH_Dfname=%s\n",MH_Dfname); |
fprintf(stderr,"MH_Dfname=%s\n",MH_Dfname); |
} |
} |
if (MH_RANK != mypower(2,MH_Mg)) { |
if (MH_RANK != mypower(2,MH_Mg)) { |
fprintf(stderr,"rk.c MH_RANK=%d is not equal to 2^MH_Mg=%d. input_data_file is broken.\n",MH_RANK,mypower(2,MH_Mg)); |
fprintf(stderr,"rk.c MH_RANK=%d is not equal to 2^MH_Mg=%d. input_data_file is broken.\n",MH_RANK,mypower(2,MH_Mg)); |
mh_usage(); |
mh_usage(); |
mh_exit(-1); |
mh_exit(-1); |
} |
} |
if (MH_Dfname != NULL) { |
if (MH_Dfname != NULL) { |
if (!MH_RawName) { |
if (!MH_RawName) { |
sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"-m%d-n%lf-x0%lf-b",MH_Mg,*MH_Ng,MH_X0g); |
sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"-m%d-n%lf-x0%lf-b",MH_Mg,*MH_Ng,MH_X0g); |
for (i=0; i<MH_Mg; i++) { |
for (i=0; i<MH_Mg; i++) { |
sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"%lf,",MH_Beta[i]); |
sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"%lf,",MH_Beta[i]); |
} |
} |
sprintf(&(MH_Dfname[strlen(MH_Dfname)]),".txt"); |
sprintf(&(MH_Dfname[strlen(MH_Dfname)]),".txt"); |
} |
} |
Line 66 mh_gopen_file() { |
|
Line 66 mh_gopen_file() { |
|
} |
} |
}else return(1); |
}else return(1); |
if (MH_byFile) { |
if (MH_byFile) { |
sprintf(fname,"%s-gp.txt",MH_Gfname); |
sprintf(fname,"%s-gp.txt",MH_Gfname); |
fp = fopen(fname,"w"); |
fp = fopen(fname,"w"); |
fprintf(fp,"set xrange [0:20]\n"); |
fprintf(fp,"set xrange [0:20]\n"); |
fprintf(fp,"set yrange [0:1.2]\n"); |
fprintf(fp,"set yrange [0:1.2]\n"); |
fprintf(fp,"plot \"%s\" title \"by hg\" with lines\n",MH_Gfname); |
fprintf(fp,"plot \"%s\" title \"by hg\" with lines\n",MH_Gfname); |
fclose(fp); |
fclose(fp); |
} |
} |
return(0); |
return(0); |
} |
} |
/* |
/* |
Runge-Kutta |
Runge-Kutta |
y_{n+1} = y_n + 1/6 k_1 + 1/3 k_2 + 1/3 k_3 + 1/6 k_4 |
y_{n+1} = y_n + 1/6 k_1 + 1/3 k_2 + 1/3 k_3 + 1/6 k_4 |
k_1 = h f(x_n, y_n) |
k_1 = h f(x_n, y_n) |
k_2 = h f(x_n + 1/2 h, y_n + 1/2 k_1) |
k_2 = h f(x_n + 1/2 h, y_n + 1/2 k_1) |
k_3 = h f(x_n + 1/2 h, y_n + 1/2 k_2) |
k_3 = h f(x_n + 1/2 h, y_n + 1/2 k_2) |
k_4 = h f(x_n + h, y_n + k_3) |
k_4 = h f(x_n + h, y_n + k_3) |
*/ |
*/ |
|
|
/* ベクトル値関数 rk is defined in tmp-code-?.c, |
/* ベクトル値関数 rk is defined in tmp-code-?.c, |
f(x, y) (x : scalar, y : vector) */ |
f(x, y) (x : scalar, y : vector) */ |
|
|
static void show_v(double x,double *v, int n) |
static void show_v(double x,double *v, int n) |
{ |
{ |
Line 109 struct MH_RESULT mh_rkmain(double x0,double y0[],doubl |
|
Line 109 struct MH_RESULT mh_rkmain(double x0,double y0[],doubl |
|
struct MH_RESULT result; |
struct MH_RESULT result; |
extern int MH_deallocate; |
extern int MH_deallocate; |
/* |
/* |
double y[MH_RANK]; |
double y[MH_RANK]; |
double k1[MH_RANK], k2[MH_RANK], k3[MH_RANK], k4[MH_RANK]; |
double k1[MH_RANK], k2[MH_RANK], k3[MH_RANK], k4[MH_RANK]; |
double temp[MH_RANK]; |
double temp[MH_RANK]; |
double ty[MH_RANK]; |
double ty[MH_RANK]; |
*/ |
*/ |
static double *y,*k1,*k2,*k3,*k4,*temp,*ty; |
static double *y,*k1,*k2,*k3,*k4,*temp,*ty; |
if (MH_deallocate && initialized) { |
if (MH_deallocate && initialized) { |
if (y) mh_free(y); |
if (y) mh_free(y); |
if (k1) mh_free(k1); |
if (k1) mh_free(k1); |
if (k2) mh_free(k2); |
if (k2) mh_free(k2); |
if (k3) mh_free(k3); |
if (k3) mh_free(k3); |
if (k4) mh_free(k4); |
if (k4) mh_free(k4); |
if (temp) mh_free(temp); |
if (temp) mh_free(temp); |
if (ty) mh_free(ty); |
if (ty) mh_free(ty); |
y = k1 = k2 = k3 = k4 = temp = ty = NULL; |
y = k1 = k2 = k3 = k4 = temp = ty = NULL; |
initialized=0; |
initialized=0; |
return(result); |
return(result); |
} |
} |
if (!initialized) { |
if (!initialized) { |
y = (double *)mh_malloc(sizeof(double)*MH_RANK); |
y = (double *)mh_malloc(sizeof(double)*MH_RANK); |
Line 135 struct MH_RESULT mh_rkmain(double x0,double y0[],doubl |
|
Line 135 struct MH_RESULT mh_rkmain(double x0,double y0[],doubl |
|
k4 = (double *)mh_malloc(sizeof(double)*MH_RANK); |
k4 = (double *)mh_malloc(sizeof(double)*MH_RANK); |
temp = (double *)mh_malloc(sizeof(double)*MH_RANK); |
temp = (double *)mh_malloc(sizeof(double)*MH_RANK); |
ty = (double *)mh_malloc(sizeof(double)*MH_RANK); |
ty = (double *)mh_malloc(sizeof(double)*MH_RANK); |
initialized=1; |
initialized=1; |
} |
} |
h = MH_Hg; |
h = MH_Hg; |
for (i = 0; i < MH_RANK; i++) |
for (i = 0; i < MH_RANK; i++) |
y[i] = y0[i]; |
y[i] = y0[i]; |
for (x = x0; (h>0?(x<xn):(x>xn)); x += h) { |
for (x = x0; (h>0?(x<xn):(x>xn)); x += h) { |
if (Df) show_v(x,y, MH_RANK); |
if (Df) show_v(x,y, MH_RANK); |
if (Gf) { |
if (Gf) { |
sprintf(swork,"%lf %le\n",x,y[0]); |
sprintf(swork,"%lf %le\n",x,y[0]); |
mh_fputs(swork,Gf); |
mh_fputs(swork,Gf); |
} |
} |
/* Output 95% point */ |
/* Output 95% point */ |
if (MH_P95) { |
if (MH_P95) { |
if ((MH_P95==1) && (y[0] >= 0.95)) { |
if ((MH_P95==1) && (y[0] >= 0.95)) { |
printf("x=%le, y[0]=%lf\n",x,y[0]); |
printf("x=%le, y[0]=%lf\n",x,y[0]); |
MH_P95=2; |
MH_P95=2; |
}else if ((MH_P95==2) && (y[0] >=0.9500000001)) { |
}else if ((MH_P95==2) && (y[0] >=0.9500000001)) { |
printf("x=%le, y[0]=%lf\n",x,y[0]); |
printf("x=%le, y[0]=%lf\n",x,y[0]); |
MH_P95=0; |
MH_P95=0; |
} |
} |
} |
} |
mh_rf(x, y, MH_RANK, temp, MH_RANK); |
mh_rf(x, y, MH_RANK, temp, MH_RANK); |
for (i = 0; i < MH_RANK; i++) |
for (i = 0; i < MH_RANK; i++) |
k1[i] = h * temp[i]; |
k1[i] = h * temp[i]; |