马文淦计算物理习题.pdf

计算物理作业 李程 PB09009006 第二章 1、采用线性同余法产生伪随机数。a = 5, c = 1, m = 16, x0 = 1。记录前 20 个数，其周期是多 少？ 前 20 个数为 1, 6, 15, 12, 13, 2, 11, 8, 9, 14, 7, 4, 5, 10, 3, 0, 1, 6, 15, 12 周期为 16 /本程序在 Visual C+ 6.0 下编译通过 #include void main() FILE *fp; int i; static int a = 5; static int c = 1; static int m = 16; int x20; x0 = 1; for(i = 1; i 0; N-) xi0 = rand(); xi = (double)(xi0 % 32767)/32767.; eta = -log(xi)/lambda; fprintf(fp,“%.15fn“,eta); fclose(fp); 6、 轻子的寿命为 3.510(-13)s，写出 N 个 轻子在实验室系中以速度 v 运动飞行距离的 抽样程序框图和程序 /本程序在 Visual C+ 6.0 下编译通过 #include #include #include #include double t = 3.5e-13; double v = 100000.0; void main() int xi0; double xi; double eta; double N; double temp; FILE *fp; printf(“请输入 N：”); scanf(“%lf”, fp = fopen(“6.txt“,“w“); srand(time(0); for (;N 0; N-) xi0 = rand(); xi = (double)(xi0 % 32767) / 32767.0; temp = v * t; eta = - temp * log(xi); fprintf(fp,“%.20fn“,eta); fclose(fp); 7、各项同性分布的角度 ， 的抽样程序 /本程序在 Visual C+ 6.0 下编译通过 #include #include #include #include void main() int xi10, xi20, xi30; double xi1, xi2, xi3; double theta, phi; double N; FILE *fp; srand(time(0); fp = fopen(“7.txt“,“w“); for (N = pow(10,6); N 0; N -) xi10 = rand(); xi1 = (double)(xi10 % 32767) / 32767.0; theta = asin(xi1); xi20 = rand(); xi2 = (double)(xi20 % 32767) / 32767.0; if (xi2 0; N-) xi10 = rand(); xi1 = (double)(xi10 % 32767) / 32767.0; x = pow(xi1,-1/n); xi20 = rand(); xi2 = (double)(xi20 % 32767) / 32767.0; y = - 1.0 / x * log(xi2); fprintf(fp,“%14.10f%14.10fn“,x,y); fclose(fp); 9、证明 ( ) ( ) ( ) 证明： 对于 ( ) 有 ( ) 对于 有 ( ) ( ) ( ) 11、正则高斯抽样 ( ) ( ) ( ) =3.2,=1.657 /本程序在 Visual C+ 6.0 下编译通过 #include #include #include void main() int xi10, xi20, xi30; double xi1, xi2, xi3; double x; double sigma, mu; int n; double N; double t1,t2; FILE *fp; printf(“请输入 sigma：“); scanf(“%lf“, printf(“请输入 mu：“); scanf(“%lf“, n = 6; N = pow(10,n); fp = fopen(“11.txt“,“w“); srand(time(0); for (; N 0; N -) xi10 = rand(); xi1 = (double)(xi10 % 32767) / 32767.0; xi20 = rand(); xi2 = (double)(xi20 % 32767) / 32767.0; t1 = pow(sigma * sigma * log(xi1) + 1.0),2); t2 = -2 * sigma * sigma * log(xi2); if (t1 = 1.0) fprintf(fp_output, “%.13fn“, x2); n_accept = n_accept + 1.; sum = sum + x2; quadsum = quadsum + x2 * x2; variance = (sum * sum - quadsum) / n_accept; else xi = uniform(0,1); if (xi 1e-4 fclose(fp_output); 6. #include #include #include #include double g(double x) return (exp(-x*x/2.0); double f(double x) return(2.506628274 * x * x); double uniform(double a, double b) int ra0; double ra; double re; ra0 = rand(); ra = (double)(ra0 % 32767) / 32767.; re = (b - a) * ra + a; return(re); void main() FILE *fp; / FILE *rand; double n_accept, n_unaccept; double x1, x2; double eta, xi; double delta; double r; double I, sum; double n; delta = 2.7; fp = fopen(“3_6.txt“,“w“); / rand = fopen(“ran.txt“,“w“); for(n = 1; n = 1.0) n_accept = n_accept + 1.; sum = sum + f(x2); / fprintf(rand,“%.13fn“,x2); else xi = uniform(0,1); if (xi = r) n_accept = n_accept + 1.; sum = sum + f(x2); / fprintf(rand,“%.13fn“,x2); else n_unaccept = n_unaccept + 1.; x2 = x1; while(n_accept pow(10,n); I = sum / (n_accept-1); fprintf(fp,“%20.13fn“,I); fclose(fp); / fclose(rand); 第七章 5. (* :Package: *) BeginPackage“legendre“; Begin“Private“ Legendren_,x_:=ModuleP,Ifn=0,P0=1,Ifn = 1,P1=x, Pn = Simplify(2n -1)x Legendren-1,x-(n-1)Legendren-2,x)/n; End EndPackag