C による数値計算法入門 ( 第 2 版 ) 新装版 サンプルページ この本の定価 判型などは, 以下の URL からご覧いただけます. http://www.morikita.co.jp/books/mid/009383 このサンプルページの内容は, 新装版 1 刷発行時のものです.
i 2 22 2 13 ( ) 2 (1) ANSI (2) 2 (3) Web http://www.morikita.co.jp/books/mid/009383 (4) C 2 Web 2 Windows 7 CPad for Boland C++ Compiler Ver.2.31 2 ( ) 2015 9 2 10 C BASIC
ii ( ) 2 ( ) ( ) C C C (1) C ANSI (2) C BASIC C C C (3) BASIC GNUPLOT (4) C Sun Solaris 2.6 gcc (GNU C version 2.8.1) ( 5 5 )
iii C 2 ( ) 2002 9 ( 1 ) ( ) (1) (2) (3) BASIC (4) (5) ( ) N88BASIC
iv 1993 1
v 1 1 1.1 2 2 1.2 6 1 9 2 1 10 2.1 1 10 2.2 1 12 2.3 15 2.4 20 2.5 1 23 2.6 28 2.7 LU 1 31 2 40 3 43 3.1 43 3.2 48 3.3 55 3.4 57 3 59 4 61 4.1 61 4.2 2 70 4 80 5 82 5.1 82 5.2 85 5.3 87 5.4 90
vi 5 93 6 94 6.1 94 6.2 97 6.3 99 6.4 2 107 6.5 2 113 6 118 7 120 7.1 120 7.2 2 127 7 129 8 130 8.1 130 8.2 132 8.3 133 8 148 9 149 9.1 149 9.2 154 9.3 160 9 169.................................................... 171........................................................ 177.......................................................... 178 178 179........................................................ 180
1 f(x) = 0 2 1.1 1.1 O 1 O 2 O 3 O 1 1 m 3 m O 3 O 1 3 O 2 2 O 2 1.1 O 2 x [m] 3 4 3 π + 2 4 3 πx3 = 4 π(3 x)3 3 x 3 3x 2 + 9x 8 = 0 (1.1) x 3 f(x) = 0
2 1 1.1 f(x) [a, b] f(a) f(b) a b c f(c) = 0 (a < c < b) 1.1 2 1.1 (1.1) 2 (1.1) f(x) f(x) = x 3 3x 2 + 9x 8 f(x) x x 1.2 f(1) = 1 < 0 f(2) = 6 > 0 [1, 2] f(x) = 0 1 2 1.5 f(1.5) f(1.5) = 2.125 > 0 f(1) < 0 f(1.5) > 0 [1, 1.5] [1, 2] 1.2
1.1 2 3 [1, 1.5] 2 2 2 1.25 f(1.25) f(1.25) = 0.51 > 0 [1, 1.25] 3 2 2 f(x) f(1.125) = 0.24 < 0, 4 : [1.125, 1.25] f(1.1875) = 0.13 > 0, 5 : [1.125, 1.1875] f(1.15625) = 0.05 < 0, 6 : [1.15625, 1.1875] f(1.171875) = 0.03 > 0, 7 : [1.15625, 1.171875] f(1.164063) = 0.01 < 0, 8 : [1.164063, 1.171875] f(1.167969) = 0.01 > 0, 9 : [1.164063, 1.167969] f(1.166016) = 0.0006 > 0, 10 : [1.164063, 1.166016] f(1.165040) = 0.005 < 0, 11 : [1.165040, 1.166016] f(1.165528) = 0.002 < 0, 12 : [1.165528, 1.166016] f(1.165772) = 0.0008 < 0, 13 : [1.165772, 1.166016] f(1.165894) = 0.00007 < 0, 14 : [1.165894, 1.166016] f(1.165955) = 0.0003 > 0, 15 : [1.165894, 1.165955] 8 x = 1.16 x = 1.17 9 x = 1.16 x mm 15 x = 1.165 x = 1.166 O 2 1 m 16 cm 6 mm 2 f(x) = 0 f(x) f(x) a 1, b 1, (a 1 < b 1 ) [a 1, b 1 ] a 1 b 1 2 c 1 f(c 1 ) f(c 1 ) = 0 c 1 ( ) f(c 1 ) 0 1.3 [a 1, c 1 ] [c 1, b 1 ] ( )
4 1 1.3 1.4 f(a 1 ) f(c 1 ) [a 1, c 1 ] 2 [a 1, c 1 ] ( 1.3(a))f(a 1 ) f(c 1 ) f(c 1 ) f(b 1 ) [c 1, b 1 ] ( 1.3(b)) 2 [c 1, b 1 ] 2 [a 2, b 2 ] [a 2, b 2 ] 3 [a 3, b 3 ] ( 1.4) 2 1.1
1.1 2 5 1.1 1 /* ************************************************* */ 2 /* 2 nibun.c */ 3 /* ************************************************* */ 4 # include < stdio.h> 5 /* ** ** */ 6 # define FNF (x) (x*x*x - 3* x*x + 9* x - 8) 7 int main ( void ) 8 { double a, b, c; 9 int k; 10 char z, zz; 11 while ( 1 ) { 12 printf ("f(a)*f(b)<0 a, b "); 13 printf (" \n\n"); 14 printf (" 1 [ a,b] a="); 15 scanf ("%lf%c",&a,& zz ); 16 printf (" 1 [ a,b] b="); 17 scanf ("%lf%c",&b,& zz ); 18 printf ("\ n (y/n)"); 19 scanf ("%c%c",&z,& zz ); 20 if(z == n ) continue ; 21 if ((z == y )&&( a < b )&&( FNF (a ) * FNF (b ) < 0)) 22 break ; 23 else { 24 printf ("\na >b f(a)*f(b ) >=0 \ n"); 25 printf (" \n"); 26 printf (" \n"); 27 scanf ("%c",&z ); continue ; 28 } 29 } 30 k = 0; 31 printf (" A B B-A\n"); 32 /* ** * * */ 33 while ( b - a >= 0.000001) { 34 k = k + 1; 35 printf ("%4d %8.5 lf %8.5 lf %8.5 lf\n",k,a,b,b-a); 36 /* a,b a,b */ 37 c = ( a + b ) / 2.0; 38 if ( FNF (a ) * FNF (c ) > 0 ) a = c; 39 else b = c; 40 if (( k % 10) == 0) { 41 printf ("\ n \n"); 42 printf (" (y/n ): "); 43 scanf ("%c%c",&z,& zz ); 44 if(z == n ) { 45 printf ("\ n \n" ); break ; 46 } else if( z == y ) 47 printf (" A B B-A\n"); 48 else { z = n ; break ; }
6 1 49 } 50 } 51 if(z!= n ) { 52 printf ("\n %3 d \ n",k); 53 printf ("\ n = %10.6 lf\n", (a+b )/2.0); 54 } 55 return 0; 56 } 1.2 (1.1) x 3 3x 2 + 9x 8 = 0 (1.1 ) f(x) [1, 2] f (x) = 3x 2 6x + 9 = 3(x 1) 2 + 6, f (x) = 6x 6 (1, 2) f (x) > 0 f (x) > 0 y = f(x) [1, 2] y = f(x) 1.5 1.5 (2, 6) x x 1 x 1 1 (x 1, f(x 1 )) x x 2 (x 2, f(x 2 )) x x 3 x 1, x 2, x 3,, x n,
40 2 1 2 0 0 0 L = 4 10 0 0 2 12 3.2 0, 1 9 6.4 6 1 4 1 1.5 1 [ U b ] = 0 1 0.6 0.5 0.4 0 0 1 1.25 0.5625 0 0 0 1 0.33333 b = t [1, 0.4, 0.5625, 0.33333] U b Ux = b u = 0.33333, z = 0.97917, y = 0.02084, x = 0.60419 LU 2 2.1 (2.4) (2.5) 1 1 2.2 (1) 2x 4y + 6z = 1 (2) 2x + 8y + 2z 3w = 2 x + 7y 8z = 0 4x + 6y 2z w = 1 x + y 2z = 3 2x 4y 2z w = 3 x 5y + 2z + w = 2 (3) 2x + 7y z + 5u 3w = 6 x 4y + 2z u + 6w = 1 3x + y 9z 2u + w = 2 10x 2y 5z + 8u 7w = 4 4x + 3y + 12z 4u 2w = 10 2.3
(1) x + 3y 2z = 2 3x 2y + z = 0 2x + y 3z = 1 (2) 3y + 2z + u = 7 3x + 2y 3z = 1 x + 2y 3z + 2u = 3 3x + 4y + z + 2u = 9 2 41 2.4 ( 2.2) 2.5 (1) 1 0 1 (2) 1 2 0 2 1 0 1 3 0 3 1 2 0 1 1 (3) 1 2 1 2 5 2 3 5 8 2.6 2.4 2.7 2.5 2.6 (1) x + y z = 0 (2) 2x 2y z = 0 (3) 5x + 2y 6z = 0 x 3y 2z = 0 4x + y 5z = 0 6x + 2y + 3z = 0 x + 2y 2z = 0 2x 2y + 2z = 0 x + 6y 6z = 0 2.8 2.8 2.9 { (1) 2x + y 4z + 5u = 1 (2) 2x 6y + 4z = 0 3x + y 2z + u = 3 x + 3y 2z = 0 3x 9y + 6z = 0 (3) x + y z = 2 (4) 2x 2y z = 12 5x + 2y 6z = 5 x 3y 2z = 19 (5) (7) 4x + y 5z = 3 x + 2y 2z = 2 2x 2y + 2z = 3 x + 6y 6z = 1 (6) 6x + 2y + 3z = 16 x + y + z = 1 2x + y 2z = 3 x y + 2z = 0 3x y + 2z = 2 x + y 2z = 2 x 1 + 2x 2 x 4 + x 6 7x 7 = 4 x 1 + 3x 2 x 3 + 3x 4 + 2x 5 + 3x 7 = 3 3x 1 + 7x 2 + x 3 + 2x 4 + 4x 5 + 3x 6 6x 7 = 13 x 1 + 2x 2 + x 4 + 4x 5 + 3x 6 11x 7 = 20 2.9 3 P Q R 3 A B C A 1 kg P 1 kg Q 1 kg R 2 kg B 1 kg P 2 kg Q 3 kg R 2 kg C 1 kg P 1 kg Q 2 kg R 3 kg A B C 1 kg 2 2.5 3
120 f(x, y) dy dx = f(x, y), y(x 0) = y 0 y(x) 1 f(x, y) y(x) x 0, x 1, x 2,, x n y(x 0 ), y(x 1 ), y(x 2 ),, y(x n ) 7.1 1 dy dx = f(x, y), y(x 0) = y 0 (7.1) x 1 = x 0 + h y 1 y(x 0 + h) y(x 0 + h) = y(x 0 ) + y (x 0 )h + 1 2! y (x 0 )h 2 + (7.2) (7.1) y (x 0 ) = f(x 0, y 0 ) f(x 0, y 0 ) y(x 1 ) y 1 (7.2) 2 y 1 = y 0 + f(x 0, y 0 )h x 1 h x 2 y(x 2 ) = y(x 1 + h) y 2 x 1 y 1 y 2 = y 1 + f(x 1, y 1 )h x x n = x 0 + nh y(x n ) y n (Euler)
7.1 121 7.1 dy dx = f(x, y), y(x 0) = y 0 h x n = x 0 + nh y(x n ) y n y n = y n 1 + f(x n 1, y n 1 )h, (n = 1, 2, ) 7.1 7.1 (7.2) h 1 h x n = x 0 + nh y(x n ) y n 7.1 y = y 12x + 3, y(0) = 1 h = 0.1 y 1, y 2, y 3, y 4, y 5, y 6, y 7, y 8, y 9, y 10 0 x 1
122 7 x j y j f(x, y) hf(x, y) x j y j f(x, y) hf(x, y) 0.0 1.0000 4.0000 0.4000 0.5 2.1159 0.8841 0.0884 0.1 1.4000 3.2000 0.3200 0.6 2.0275 2.1725 0.2173 0.2 1.7200 2.3200 0.2320 0.7 1.8102 3.5898 0.3590 0.3 1.9520 1.3520 0.1352 0.8 1.4512 5.1488 0.5149 0.4 2.0872 0.2872 0.0287 0.9 0.9363 6.8637 0.6864 1.0 0.2499 8.7501 0.8750 7.2 1 y = 12x 8e x + 9 x = 0.1 0.5 1 y(0.1) = 1.3586, y(0.5) = 1.8102, y(1.0) = 0.7463 (7.2) h 1 2 3 4 O(h n ) n h k h 0 k/h n c( ) k = O(h n ) ( c = 0 ) O(h n ) + O(h n ) = O(h n ), ho(h n ) = O(h n+1 ) 0 < n < m O(h n ) + O(h m ) = O(h n ) (7.3) y(x) C n+1 y(x) y(x 0 + h) = y(x 0 ) + y (x 0 )h + 1 2! y (x 0 )h 2 + + 1 n! y(n) (x 0 )h n + O(h n+1 ) f(x, y) y f(a, b + h) = f(a, b) + O(h) (7.4) y(x 0 +h) h n y(x 0 +h) h n y(x 0 ) = y 0 y (x 0 ) = y 0 y (x 0 ) = y 0
7.1 123 h 2 y(x 0 + h) = y 0 + y 0h + 1 2! y 0 h 2 + O(h 3 ) (7.5) (7.5) y 0 = f(x 0, y 0 ) y 0 y 0 f(x, y) (7.5) h 2 f(x, y) y = y(x 0 + h) y(x 0 ) = hy (x 0 + θh), (0 < θ < 1) y (x 0 + θh) θ = 0 θ = 1 θ = 0 y hy (x 0 ) θ = 1 y hy (x 0 + h) y h 2 αhy (x 0 ) + βhy (x 0 + h) α β y h 2 αhy (x 0 ) + βhy (x 0 + h) = αhy (x 0 ) + βh{y (x 0 ) + y (x 0 )h + O(h 2 )} = (α + β)hy 0 + βh 2 y 0 + O(h 3 ) (7.5) α + β = 1, β = 1 2, α = 1 2, β = 1 2 y = 1 2 hy (x 0 ) + 1 2 hy (x 0 + h) + O(h 3 ) k 1 = hy (x 0 ) = hf(x 0, y 0 ) (7.6) (7.5) (7.4) hy (x 0 + h) = hf(x 0 + h, y(x 0 + h)) = hf(x 0 + h, y(x 0 ) + y (x 0 )h + O(h 2 )) = hf(x 0 + h, y 0 + k 1 + O(h 2 )) = hf(x 0 + h, y 0 + k 1 ) + O(h 3 ) k 2 = hf(x 0 + h, y 0 + k 1 ) (7.7)
124 7 y = 1 2 k 1 + 1 2 k 2 + O(h 3 ) k = 1 2 (k 1 + k 2 ), y 1 = y 0 + k (7.8) y(x 0 + h) = y 1 + O(h 3 ) y 1 h 2 y 2 x 1 y ( ) y 1 (7.6) (7.7) (7.8) k 1 = hf(x 1, y 1 ), k 2 = hf(x 1 + h, y 1 + k 1 ), k = 1 2 (k 1 + k 2 ) y 2 = y 1 + k y 3 y 4 (Runge-Kutta) 2 7.2 2 dy dx = f(x, y), y(x 0) = y 0 h x n = x 0 + nh y n y n+1 y n+1 = y n + k, (n = 0, 1, 2, ) k k 1 = hf(x n, y n ), k 2 = hf(x n + h, y n + k 1 ), k = 1 2 (k 1 + k 2 ) 2 7.1 x j y j k j k 0.0 1.0 0.4 0.1 1.4 0.32 0.36 0.1 1.36 0.316 0.2 1.676 0.2276 0.2718 0.2 1.6318 x j y j 0.3 1.80614 0.4 1.87279 0.5 1.82043 0.6 1.63657 x j y j 0.7 1.30741 0.8 0.81769 0.9 0.15055 1.0 0.71265
7.1 125 7.2 7.1 7.1 1 /* ************************************************* */ 2 /* 2 rungekt2.c */ 3 /* ************************************************* */ 4 # include <stdio.h> 5 double fnf ( double x, double y) 6 { return ( y - 12.0 * x + 3.0); } 7 int main ( void ) 8 { int i; 9 double x, y, h, k1, k2, k; 10 char zz; 11 printf (" 2 \ n\n"); 12 printf ("dy/dx = y - 12.0 * x + 3.0 "); 13 printf ("\n\ n \n"); 14 scanf ("%c",&zz ); 15 printf (" X Y\n"); 16 x = 0.0; y = 1.0; h = 0.1; 17 printf (" %10.6 lf %10.6 lf\n",x,y); 18 for (i =1; i <=20; i ++) { 19 k1 = h * fnf (x,y); 20 k2 = h * fnf (x+h,y+k1 ); 21 k = ( k1 + k2 ) / 2.0; 22 x = x + h; 23 y = y + k; 24 printf (" %10.6 lf %10.6 lf\n",x,y); 25 } 26 return 0; 27 } 2 h 4 4 ( [17]) 7.3 4 dy dx = f(x, y), y(x 0) = y 0 h x n = x 0 + nh y n y n+1 = y n + k, (n = 0, 1, 2, ) k
171 (2-1) (2-1) 2 1 1 1.1 (1) 0.3376 1.3075 (2) 0.7390 (3) 0.5671 1.2 0.53727 1.31597 1.3 1.4902 1.4 21.3 cm 1.5 t = 2.4587 2 x 1 = b 1 a 21 x 1 + x 2 = b 2 x 1 = b 1 j 1 2.1 a 31 x 1 + a 32 x 2 + x 3 = b 3,. x j = b j a jk x k. k=1 a n1 x 1 + a n2 x 2 + + a n,n 1 x n 1 + x n = b n (j = 2, 3,, n) 2.2 (1) x = 1.8 y = 1 z = 1.1 (2) x = 0.166667 y = 0.166667 z = 0.541667 w = 1.583333 (3) x = 0.736418 y = 0.014627 z = 0.349701 u = 1.837313 w = 0.721940 2.3 (1) x = 0.375 y = 0.625 z = 0.125 (2) x = 2 y = 1 z = 3 u = 4 2.4 (2-1) 2 1 1 2.5 (1) 4 1 2 1 1 1 3 (2) 2.6 (2-2) 4 2.7 (1) x = t 1 (2) x = t (t ) 3 2 0 1 1 0 1 1 1 4 (3) 1 0 3 (3) x = t 1 50 11 9 22 5 4 5 1 1 2.8 s t x 1
172 4 2 4 3 2 (1) x = 9 0 + s 8 1 + t 13 0 (2) x = s 1 + t 0 0 1 0 0 1 (3) x = 1 1 4 0.5 1 5 + s 1 (4) x = 6.5 + t 3 (5) 3 0 3 0 4 84 14 12 52 36 6 5 21.5 1 3 0 0 3.5 (6) x = s 1 (7) x = 8 + r 2 + s 1 + t 2 1 0 1 0 0 0 0 1 0 0 0 0 1 2.9 x = 11 3, y = 1 3, z = 5 79 kg, 3 6 2 0 0 1 3 4 2.10 (1) 1 4 0 0 1 2 = 40 (2) (3) 3 5 5 2 0 0 3 4 0 3 2 1 0 0 1 1 2 0 0 1 3 = 8 1 5 1 0 0 1 2 0 0 1 3 5 1 4 0 0 1 2 = 8 0 0 1 2.11 2.12 3 0 0 0 1 2 1 3 2.13 (1) 1 2 0 0 0 1 3 4 1 4 4 0 0 0 1 2.5 = 96 0 2 4 4 0 0 0 1 3 0 0 0 1 0.5 2 1.6 (2) 1 1 0 0 0 1 0 4 0 1.5 2 0 0 0 1 3.5 = 3 2 3 2.2 0.5 0 0 0 1
178 1.1 2 5 1.2 8 2.1 14 2.2 18 2.3 LU 37 3.1 47 3.2 54 4.1 (1 ) 67 4.2 2 78 6.1 96 6.2 106 6.3 108 6.4 1 109 6.5 2 111 6.6 2 116 7.1 2 125 8.1 135 8.2 ( ) 141 9.1 158 9.2 () 166
179 [a 1, a 2,, a n ] 10 a 1 a 2 10. a n [a ij ] a ij i j 11 t a a 12 t A A 75 A 1 A 22 A A 37 n a j a 1 a 2 a n 45 j=1 L k (x) n j=1 j k x x j x k x j 45 f[x 0, x 1 ] 1 49 f[x 0, x 1, x 2 ] 2 49 f[x 0, x 1,, x n ] n 50 C n n 51 ( ) 55 2 2 56 n n 56 nc j 58 ( n j ) 58 T n (x) n 82 ζ 0, ζ 1, ζ 2,, ζ n T n+1 (x) 87 P n (x) n 90 H n (k) 104 O(h n ) h n ( ) 122 x x 149
180 1... 56 1... 64 1... 131 1... 131 2... 56 2... 130 2... 127 2... 65 2... 131 2... 111 2... 113 2... 2 3... 56 C n... 51 DE... 111 h n... 122 LU... 31 n... 56 n 1... 11... 12 1... 12... 120 138... 96... 122... 20... 18... 101... 12... 39... 21... 13... 47... 133... 130... 12... 61... 138... 11... 132... 132... 151... 149... 149... 151 2... 70... 85... 52... 51... 55... 56 1... 12... 36... 133... 97... 96... 61... 75... 61... 132... 18... 131 1... 49 n... 50... 95... 160
181... 131... 83... 87... 87... 88... 2... 133... 160... 120... 50... 6... 7... 17... 154... 154 x... 149... 130... 144... 131... 46... 46... 43... 103... 160... 44... 45... 43... 144... 90 2... 124 4... 125 4 128