x(t) + t f(t, x) = x(t) + x (t) t x t Tayler x(t + t) = x(t) + x (t) t + 1 2! x (t) t ! x (t) t 3 + (15) Eular x t Teyler 1 Eular 2 Runge-Kutta

Size: px

Start display at page:

Download "x(t) + t f(t, x) = x(t) + x (t) t x t Tayler x(t + t) = x(t) + x (t) t + 1 2! x (t) t ! x (t) t 3 + (15) Eular x t Teyler 1 Eular 2 Runge-Kutta"

まともえいさか
6 years ago
Views:

1 6 Runge-KuttaEular Runge-Kutta Runge-Kutta A( ) f(t, x) dx dt = lim x(t + t) x(t) t 0 t = f(t, x) (14) t x x(t) t + dt x x(t + dt) Euler 7 t 1 f(t, x(t)) x(t) + f(t + dt, x(t + dt))dt t + dt x(t + dt) (14) x(t + t) = x x(t+dt) x(t) t t+dt 7: 1

2 x(t) + t f(t, x) = x(t) + x (t) t x t Tayler x(t + t) = x(t) + x (t) t + 1 2! x (t) t ! x (t) t 3 + (15) Eular x t Teyler 1 Eular 2 Runge-Kutta 8 Euler t+dt x k 1 t t+dt t+dt/2, x(t)+ k 1 /2 x f(t + dt/2, x(t) + k 1 /2) x(t + dt) x x(t+dt) x(t) k 1 t t+dt/2 t+dt 8: 2 Runge-Kutta k 1 = f(t, x(t))dt (16) x(t + dt) = x(t) + f(t + dt/2, x(t) + k 1 /2)dt (17) 2 Runge-Kutta Taylor dt 2 Runge-Kutta 9 4 Runge-Kutta (16) (17) 1. (t, x) dt x k 1 2. k 1 t t + dt f(t + dt, x + k 1 /2) x(t) dt x k 2 3. k 2 t t + dt f(t + dt, x + k 2 /2) x(t) dt x k 3 2

3 4. k 3 t t + dt f(t + dt, x + k 3 /2) x(t) dt x k 4 x k 1 k 2 k 3 k 4 x(dt) k 1 = f(t, x(t))dt (18) k 2 = f(t + dt/2, x(t) + k 1 /2)dt (19) k 3 = f(t + dt/2, x(t) + k 2 /2)dt (20) k 4 = f(t + dt, x(t) + k 3 )dt (21) x(t + dt) = x(t) (k 1 + 2k 2 + 2k 3 + k 4 ) (22) Eular t+dt (1/6)t+dt (1/2) t + dt (5/6) 3 4 Runge-Kutta Taylor 4 dt 5 Runge-Kutta 4 Runge-Kutta x x dt/6 dt/3 dt/3 dt/6 x(t+dt) x(t) k 1 k 2 k 3 k 4 step1 step2 step3 step4 t t+dt/2 t+dt t t+dt/2 t+dt 9: 4 Runge-Kutta dy/dx = y Runge-Kutta Runge-Kutta S_RungK 6 Runge-Kutta // === Runge-Kutta === #include <stdio.h> #include <math.h> #define EPS 1.0e-8 #define N 50 // 3

4 #define X0 0.0 // X #define XN 5.0 // X #define Y0 1.0 // Y // main main // // void S_RungK( double, double *, double, int, double(*)() ); double FUNC( double, double ); // void main(){ FILE *fout; int i; double x, h, y[101]; fout=fopen("rk.csv","w"); h=(xn-x0)/(double)n; x=x0; y[0]=y0; printf(" X0=%f XN=%f Y0=%f h=%f \n",x,xn,y[0],h); S_RungK(x,y,h,N,*FUNC); for(i=0; i<=n; i++) { printf(" i=%f, y=%f \n", h*i, y[i]); fprintf(fout,"%f, %f\n", h*i, y[i]); fclose(fout); /** Function to calculate Runge-Kutta method ** // ; y[0] = y // h = x // N = N // x = x // FUNC = dy/dx // ; y[n+1=n+1 void S_RungK( double x, double y[], double h, int n, double (*FUNC)()){ int i; double k1,k2,k3,k4; for(i=0; i<n; i++){ k1=h*(*func)(x, y[i]); // step1 k2=h*(*func)(x+0.5*h, y[i]+0.5*k1); // step2 k3=h*(*func)(x+0.5*h, y[i]+0.5*k2); // step3 k4=h*(*func)(x+h,y[i]+k3); // step4 x=x+h; y[i+1]=y[i]+(k1+2.0*(k2+k3)+k4)/6.0; 4

5 // dy/dx // dy/dx=y double FUNC( double x, double y ){ return ( y ); y : dy/dx = y 4 Runge-Kutta x Runge-Kutta qx, qy, px, py RungeKutta() 7 Runge-Kutta // === 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 5

6 7 2 m [kg] q x, q y q y dq 2 x dt = 0, dqy 2 = g (23) 2 dt2 2 x y p x [kg m/s] p y [kg m/s] dq 2 x dt = dp x 2 dt = 0, dq2 y dt = dp y = g, 2 dt (24) dq x dt = p x m, dq y dt = p y m, (25) dt p x p y 1 p x [m]p y [m] 1 Eular 8 Eular // Eular #include <stdio.h> #define G // [m/s/s] #define M 1.0 // [kg] #define N 256 // #define T0 0.0 // [s] #define TN 5.0 // [s] #define QX0 0.0 // X [m] #define QY0 0.0 // Y [m] #define PX0 1.0 // X [kgm/s] #define PY0 2.0 // Y [kgm/s] //---- main function int main(void) { FILE *fout; int i; double qx, qy, px, py, dt; fout = fopen("kekka.csv","w"); // printf(" T0=%4.1f QX0%4.1f QY0=%4.1f PX0=%4.1f PY0=%4.1f N=%5d\n",T0,QX0,QY0,PX0,PY0,N); 6

7 dt = (TN-T0)/(double)N; qx = QX0; qy=qy0; px=px0; py=py0; for(i=0; i<n; i++) { // qy printf("qx=%lf yx=%lf\n",qx,qy); // qx,qy fprintf(fout,"%lf, %lf\n",qx,qy); // qx,qy qx += px/m*dt; // x () qy += py/m*dt; // y () px += M*0*dT; // x () py += -M*G*dT; // x () if(qy<0.0) break; // y 0 fclose(fout); return(0); // 垂直方向 y 水平方向 x 11: () Eular Runge-Kutta

8 (q x, q y ) (p x, p y ) p x = m dq x dt, p y = m dq y dt dp x dt = GMm cos ϕ = GMmq x, q 2 q 3 dp y dt = GMm q 2 (26) cos ϕ = GMmq y q 3, (27) GMm r = mrω 2, r = 1, T = 1, ω = 2π T = 2π (28) G M GM 4π 2 1 dt Runge-Kutta Runge-Kutta q x q y p x p y RungeKutta() 9 Runge-Kutta // Runge-Kutta #include <stdio.h> #include <math.h> // ---- physical setting #define M_PI #define GM 4*M_PI*M_PI #define M 1.0 // #define dt 1.0/256 // #define QX0 1.0 // x #define QY0 0.0 // y #define PX0 0.0 // x #define PY0 sqrt(gm)*m // y // ---- main function int main(void) { FILE *fout; double qx, qy; double px, py; 8

9 double T; double qxk1, qyk1, pxk1, pyk1; double qxk2, qyk2, pxk2, pyk2; double qxk3, qyk3, pxk3, pyk3; double qxk4, qyk4, pxk4, pyk4; double tqx, tqy, tpx, tpy; double q, qi3; fout = fopen("orbit.csv","w"); // qx = QX0; qy = QY0; px = PX0; py = PY0; for( T = 0.0; T < 10.00; T += dt ){ // (qx,qy) (px,py) k1 tqx = qx; tqy = qy; tpx = px; tpy = py; q = hypot( tqx, tqy ); qi3 = 1.0/(q*q*q); qxk1 = tpx / M * dt; qyk1 = tpy / M * dt; pxk1 = -GM * M * tqx * qi3 * dt; pyk1 = -GM * M * tqy * qi3 * dt; // (qx,qy) (px,py) k2 tqx = qx+0.5*qxk1; tqy = qy+0.5*qyk1; tpx = px+0.5*pxk1; tpy = py+0.5*pyk1; q = hypot( tqx, tqy ); qi3 = 1.0/(q*q*q); qxk2 = tpx / M * dt; qyk2 = tpy / M * dt; pxk2 = -GM * M * tqx * qi3 * dt; pyk2 = -GM * M * tqy * qi3 * dt; // (qx,qy) (px,py) k3 tqx = qx+0.5*qxk2; tqy = qy+0.5*qyk2; tpx = px+0.5*pxk2; tpy = py+0.5*pyk2; q = hypot( tqx, tqy ); qi3 = 1.0/(q*q*q); qxk3 = tpx / M * dt; qyk3 = tpy / M * dt; pxk3 = -GM * M * tqx * qi3 * dt; pyk3 = -GM * M * tqy * qi3 * dt; // (qx,qy) (px,py) k4 tqx = qx+qxk3; tqy = qy+qyk3; tpx = px+pxk3; tpy = py+pyk3; q = hypot( tqx, tqy ); qi3 = 1.0/(q*q*q); qxk4 = tpx / M * dt; qyk4 = tpy / M * dt; pxk4 = -GM * M * tqx * qi3 * dt; pyk4 = -GM * M * tqy * qi3 * dt; // k1,k2,k3,k4 qx += (qxk1 + 2*qxk2 + 2*qxk3 + qxk4)*(1.0/6); qy += (qyk1 + 2*qyk2 + 2*qyk3 + qyk4)*(1.0/6); px += (pxk1 + 2*pxk2 + 2*pxk3 + pxk4)*(1.0/6); py += (pyk1 + 2*pyk2 + 2*pyk3 + pyk4)*(1.0/6); printf("qx=%f qy=%f px=%f py=%f \n",qx,qy,px,py); 9

10 fprintf(fout,"%f, %f\n",qx,qy); // qx,qy fclose(fout); return(0); // 12 Runge-Kutta Euler 10 Euler (?) y : Keplar Runge-Kutta x (1) (2) (3) t x t x 10

11 u(x, t) u(x, t) t = κ 2 u(x, t) x 2 (29) κ 9.1 SI CGS T X /T τx/x ξ 1 κ = 1 u(x, t) t = 2 u(x, t) x 2 (30) 9.2 t t > 0 x 0 < x < LL 0 < x < L u(x, 0) = 1 x = 0 x = L 0u(0, t) = 0t > 0 u(x, t) tx t x u(x, t) u(j x, (n + 1) t) u(j x, n t) (31) t t 2 u(x, t) u((j + 1) x, n t) 2u(j x, n t) + u((j 1) x, n t) x 2 x 2 (32) J max = 1/ x j = 0, 1, 2,, J max n = 0, 1, 2, n u(j x, (n + 1) t) u(j x, n t) t (33) u((j + 1) x, n t) 2u(j x, n t) + u((j 1) x, n t) = x 2 (34) 11

12 t u(j x, (n + 1) t) u(j x, n t) (35) = t x {u((j + 1) x, n t) 2u(j x, n t) + u((j 1) x, n t) 2 (36) t/ x 2 r u(j t, (n + 1) x) = r u((j + 1) x, n t) +(1 2r) u(j x, n t) + r u((j 1) x, n t) (37) r = t x 2 (38) 10 // #include <stdio.h> #define JMAX 10 // #define NMAX 250 // #define NINT 25 // void main() { FILE *fout; double u[jmax+1],unew[jmax+1]; // u unew double delx,delt,r; // x t int j,n; // j n delx = 1.0/(double)JMAX; delt = 0.004; r = delt/delx*delx; // x // // r // u[0] = 0.0; for(j=1; j<jmax; j++) u[j]=1.0; u[jmax]=0.0; for(j=0; j<=jmax; j++) unew[j]=u[j]; fout = fopen("diffuse.csv","w"); // for( n=0; n<=nmax; n++) { // // NINT if( (n%nint) == 0) { for(j=0; j<jmax; j++) fprintf(fout, "%lf, ", u[j]); fprintf(fout,"%lf\n",u[jmax]); for(j=1; j<jmax; j++) unew[j] = r*u[j+1]+(1.0-r*2.0)*u[j] +r*u[j-1]; 12

13 // 1 // for(j=0; j<=jmax; j++) u[j] = unew[j]; // fclose(fout); // r < 1 2 (39) x t 13

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