x h = (b a)/n [x i, x i+1 ] = [a+i h, a+ (i + 1) h] A(x i ) A(x i ) = h 2 {f(x i) + f(x i+1 ) = h {f(a + i h) + f(a + (i + 1) h), (2) 2 a b n A(x i )

1 f(x) a b f(x)dx = n A(x i ) (1) ix [a, b] n i A(x i ) x i 1 f(x) [a, b] n h = (b a)/n y h = (b-a)/n y = f (x) h h a a+h a+2h a+(n-1)h b x 1: 1

x h = (b a)/n [x i, x i+1 ] = [a+i h, a+ (i + 1) h] A(x i ) A(x i ) = h 2 {f(x i) + f(x i+1 ) = h {f(a + i h) + f(a + (i + 1) h), (2) 2 a b n A(x i ) = = h n h 2 {f(x i) + f(x i+1 ) (3) { n 1 1 2 f(a) + f(a + i h) + 1 2 f(b) (4) 1 1 #include <stdio.h> #define X_MIN 0.0 // #define X_MAX 100.0 // #define DIV_NUM 120 // //. float function (float x) { float result; result = x; // return result; int main() { double integral; // double h; // n double x, da; int i; printf ("=[%f,%f] =%d\n", X_MIN, X_MAX, DIV_NUM); // === === h = (X_MAX - X_MIN) / DIV_NUM; // integral = 0.0; // x = X_MIN; // for (i=0; i<div_num; i++) { // A da = (function(x)+function(x+h))*h/2.0; integral += da; // 2

x += h; // === === printf ("= %lf\n", integral); printf (" ( 100 100 )=%lf\n", (X_MAX-X_MIN)*function(X_MAX)/2.0); return(0); (3) f(x + h) f(x) 2 2 2 2 y h = (b-a)/n y = f (x) h/2 h/2 a a+h a+(h/2) b a+{(2n-1)/2h x 2: f(x) [a, b] n h = (b a)/n [x i, x i+1 ] = [a + i h, a + (i + 1) h] x i + h/2 h A(x i ) A(x i ) = h f (x i + 12 ) ( h = h f a + 2i 1 ) h, (5) 2 3

a b n A(x i ) = h n f (x i + 12 ) h = h n ( f a + 2i 1 ) h 2, (6) 1 3 f(x) Simpson f(x) 2 3 (x 0, f(x 0 ))(x 1, f(x 1 ))(x 2, f(x 2 )) 3 g(x) g(x) = (x x 1)(x x 2 ) (x 0 x 1 )(x 0 x 2 ) f(x 0) + (x x 2)(x x 0 ) (x 1 x 2 )(x 1 x 0 ) f(x 1) + (x x 0)(x x 1 ) (x 2 x 0 )(x 2 x 1 ) f(x 2) (7) 2 f (x) f (x) x 0 x 2 x 0 x 1 x 2 3: () () 4

g(x) [x 0, x 2 ] x2 [ (x x1 )(x x 2 ) g(x)dx = x 0 2(x 0 x 1 )(x 0 x 2 ) f(x (x x 1 ) 3 0) 6(x 0 x 1 )(x 0 x 2 ) f(x 0) + (x x 2)(x x 0 ) 2(x 1 x 2 )(x 1 x 0 ) f(x 1) (x x 2 ) 3 6(x 1 x 2 )(x 1 x 0 ) f(x 1) + (x x 0)(x x 1 ) 2(x 2 x 0 )(x 2 x 1 ) f(x (x x 0 ) 3 2) 6(x 2 x 0 )(x 2 x 1 ) f(x 2) = (x 2 x 0 ) 2 (x 2 x 1 ) 2(x 2 x 0 )(x 2 x 1 ) f(x (x 2 x 0 ) 3 2) 6(x 2 x 0 )(x 2 x 1 ) f(x 2) + (x 0 x 1 )(x 0 x 2 ) 2 2(x 0 x 1 )(x 0 x 2 ) f(x (x 0 x 2 ) 3 0) + 6(x 0 x 1 )(x 0 x 2 ) f(x 0) (x 0 x 2 ) 3 + 6(x 1 x 2 )(x 1 x 0 ) f(x 1) (8) x 1 x 0 = x 2 x 1 = (x 2 x 0 )/h (8) x2 x 0 g(x)dx = (2h)2 h 2(2h)h f(x 2) (2h)3 6(2h)h f(x 2) ( h)( 2h)2 2( h)( 2h) f(x 0) + ( 2h)3 6( h)( 2h) f(x 0) + ( 2h)3 6( h)h f(x 1) =hf(x 2 ) 4 6 hf(x 2) + hf(x 0 ) 4 6 hf(x 0) + 8 6 hf(x 1) = h 3 [f(x 0) + 4f(x 1 ) + f(x 2 )] (9) f(x) [a, b] n A(x i ) = h 3 [f(x 0) + 4f(x 1 ) + f(x 2 )] + h 3 [f(x 2) + 4f(x 3 ) + f(x 4 )] + + h 3 [f(x 2n 2) + 4f(x 2n 1 ) + f(x 2n )] = h 3 [f(x 0) + 4f(x 1 ) + 2f(x 2 ) + 4f(x 3 ) + 2f(x 4 ) + 2f(x 2n 2 ) + 4f(x 2n 1 ) + 2f(x 2n )] = h 3 [f(x 0) + 4f(x 1 ) + f(x 3 ) + + 2f(x 2 ) + f(x 4 ) + + f(x n )] [ ] = h n n 1 f(a) + 4 fa + (2i 1)h + 2 fa + 2ih + f(b) (10) 3 2 Simpson 5 ] x2 x 0

f (x) f (x) h da 1 da 2 da n-1 da n a b x x 0 x 1 x 2 x 3 x 4 x 2n-2 x 2n-1 x 2n 4: // === Simpson ( === h = (X_MAX - X_MIN) / (2.0*DIV_NUM); // x = X_MIN; // integral = function(x); // for (; i<div_num; i++) { da = 4.0*function(x+h) + 2.0*function(x+2.0*h); integral += da; // x += 2.0*h; integral += ( 4.0*function(x+h) + function(x+2.0*h) ); integral *= h/3.0; // === Simpson ( === 4 V N x 1, x 2,, x n V f fdv V V (11) V f Monte Carlo C rand() 3 6

A f (x)dx 5: A f f A #include <stdio.h> #include <stdlib.h> // rand #define N 400 // #define D 10.0 // main() { unsigned long i; unsigned long in=0,out=0; float ratio,x,y,r,square; square=d*d; // r=d/2.0; // for(i=0; i<n; ++i) { // x x= (double)rand()/(double)rand_max*r; // y y= (double)rand()/(double)rand_max*r; // (x,y) if ((x*x+y*y)<=r*r) ++in; else ++out; printf("=%d =%d =%d x=%5.2f y=%5.2f\n", i,in,out,x,y); // (x,y) ratio=(float)in/(float)(in+out); printf("= %d = %d\n",in,out); printf("%f\n",square); printf(" %f\n",square*ratio); printf("3.14*r*r%f\n",3.14*r*r); 7

5 http://www.sm.rim.or.jp/~shishido/bibun1.html dy dx = 2x, x = 1 y = 1 (12) y = x 2 ) y dy/dx = 2x x = 1 y 1 2 x, 2 x = 1.1y x(0.1)x = 1 y x = 1 y x = 1.1 y y 2 0.1 = 0.2 y 1 1.2 (1.1)2 = 1.21 x = 1 y(x + x) = y(x) + dy dx x (13) x y x 2 1 x 1.0 6 x = 1.0 y = 1.0 h = 0.02 4 dy dx #include <stdio.h> #define XSTART 1.0 // x #define XEND 6.0 // x #define YSTART 1.0 // y #define H 0.02 // x main() { double x, y, dy; y = YSTART; x = XSTART; = 2x 8

while(x<6.0) { dy = 2.0*x*H; y += dy; x += H; printf("x=%lf y=%lf\n",x,y); return(0); x=1.020000 y=1.040000 x=1.040000 y=1.080800 x=1.060000 y=1.122400 x=1.080000 y=1.164800 x=1.100000 y=1.208000 x=5.940000 y=35.184800 x=5.960000 y=35.422400 x=5.980000 y=35.660800 x=6.000000 y=35.900000 x=6.020000 y=36.140000 8 10 5 dy dx = 2x #include <stdio.h> #include <glut.h> // OpenGL #define XSTART 1.0 // x #define XEND 6.0 // x #define YSTART 1.0 // y #define H 0.02 // x void display(void) // { double x, y, dy; glclear(gl_color_buffer_bit); // glbegin(gl_line_strip); // ( y = YSTART; x = XSTART; while(x<6.0) { // Eular dy = 2.0*x*H; y += dy; x += H; printf("x=%lf y=%lf\n",x,y); glvertex2d(x*0.2-0.9, y*0.04-0.9); // 9

glend(); // ( glbegin(gl_line_strip); // ( glvertex2d(-0.9, 0.9); glvertex2d(-0.9, -0.9); // glvertex2d( 0.9, -0.9); glend(); // ( glflush(); // 6: 1 2 2 1 1 2 6 Runge-Kutta A da dt = f(a) (14) f(a) t = 0 A A(0) t = dt A A(dT ) Euler 7 t = 0 1 f(a(0)) dt A(dT ) f(a(0))dt 10

x x(t+dt) x(t) t t+dt 7: 8 2 Runge-Kutta Euler dt A k 1 0 dt dt/2 A f(a(0)) + k 1 /2 A(dT ) k 1 = f(a(0))dt (15) A(dT ) = A(0) + f(a(0) + k 1 /2)dT (16) Taylor dt 2 Runge-Kutta 9 4 Runge- Kutta (15) (16) 4 Runge-Kutta 3 A(dT ) k 1, k 2, k 3, k 4 A(dT ) Taylor dt 4 k 1 = f(a(0))dt (17) k 2 = f(a(0) + k 1 /2)dT (18) k 3 = f(a(0) + k 2 /2)dT (19) k 4 = f(a(0) + k 3 /2)dT (20) A(dT ) = A(0) + 1 6 (k 1 + 2k 2 + 2k 3 + k 4 ) (21) Runge-Kutta 4 Runge-Kutta Runge-Kutta qx, qy, px, py RungeKutta() 6 Runge-Kutta // === Runge-Kutta === 11

x x(t+dt) t x(t) t+dt/2 t+dt k 1 8: 2 Runge-Kutta x x(t+dt) x(t) k 1 k 2 k 3 k 4 t t+dt 9: 4 Runge-Kutta void RungeKutta(double qx, double qy, double px, double py, double *qxk, double *qyk, double *pxk, double *pyk) { double q, qi3; q = hypot( qx, qy ); qi3 = 1.0/(q*q*q); *qxk = px/m*dt; *qyk = py/m*dt; *pxk = -GM*M*qx*qi3*dT; *pyk = -GM*M*qy*qi3*dT; Runge-Kutta Euler 10 Euler 12