home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / ode / rungek.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  7.1 KB  |  366 lines
  1.  
  2.  
  3. /* Runge-Kutta numerical integration
  4.  * of a system of ordinary differential equations.
  5.  *
  6.  * Reference:
  7.  * Thomas, Benku, "The Runge-Kutta Methods," BYTE 11, #4, April 1986
  8.  * This program is adapted from code supplied with that article.
  9.  */
  10.  
  11. #include "int.h"
  12. #define ORDER 7
  13. #define ERRTERM 0
  14.  
  15. #if !(ORDER-4)
  16. /*
  17. Fourth order Runge-Kutta method specified by
  18. Felberg, E., COMPUTING,    6(1970)p61-71
  19. */
  20. #define INUM 6
  21. #define HMAXL 0.5
  22. #define HLIM1 3.0
  23.  
  24. double al[INUM-1] = {
  25. 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0
  26. };
  27.  
  28. double b[] = {
  29. 1.0/4.0, /* b(1,1) */
  30.  
  31. 3.0/32.0, /* b(1,2) */
  32. 9.0/32.0, /* b(2,2) */
  33.  
  34. 1932.0/2197.0, /* b(1,3) */
  35. -7200.0/2197.0, /* b(2,3) */
  36. 7296.0/2197.0, /* b(3,3) */
  37.  
  38. 439.0/216.0, /* b(1,4) */
  39. -8.0,          /* b(2,4) */
  40. 3680.0/513.0, /* b(3,4) */
  41. -845.0/4104.0, /* b(4,4) */
  42.  
  43. -8.0/27.0,      /* b(1,5) */
  44. 2.0,             /* b(2,5) */
  45. -3544.0/2565.0, /* b(3,5) */
  46. 1859.0/4104.0, /* b(4,5) */
  47. -11.0/40.0,    /* b(5,5) */
  48. };
  49.  
  50. #if ERRTERM
  51. double er[INUM] = {
  52. 1.0/360.0, 0.0, -384.0/12825.0, -41743.0/1429560.0, 1.0/50.0, 2.0/55.0
  53. };
  54. #endif
  55.  
  56. double a[INUM] = {
  57. 16.0/135.0, 0.0, 6656.0/12825.0, 28561.0/56430.0, -9.0/50.0, 2.0/55.0
  58. };
  59. #endif
  60.  
  61.  
  62. #if !(ORDER-5)
  63. /*
  64. Fifth order Runge-Kutta method specified by
  65. Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 5)
  66. */
  67. #define INUM 8
  68. #define HMAXL 1.0
  69. #define HLIM1 4.0
  70. double al[INUM-1] = {
  71. 1.0/18.0, 1.0/6.0, 2.0/9.0, 2.0/3.0, 1.0, 8.0/9.0, 1.0
  72. };
  73.  
  74. double b[] = {
  75. 1.0/18.0, /* b(1,1) */
  76.  
  77. -1.0/12.0, /* b(1,2) */
  78. 1.0/4.0, /* b(2,2) */
  79.  
  80. -2.0/81.0, /* b(1,3) */
  81. 4.0/27.0, /* b(2,3) */
  82. 8.0/81.0, /* b(3,3) */
  83.  
  84. 40.0/33.0, /* b(1,4) */
  85. -4.0/11.0, /* b(2,4) */
  86. -56.0/11.0, /* b(3,4) */
  87. 54.0/11.0, /* b(4,4) */
  88.  
  89. -369.0/73.0, /* b(1,5) */
  90. 72.0/73.0, /* b(2,5) */
  91. 5380.0/219.0, /* b(3,5) */
  92. -12285.0/584.0, /* b(4,5) */
  93. 2695.0/1752.0, /* b(5,5) */
  94.  
  95. -8716.0/891.0, /* b(1,6) */
  96. 656.0/297.0, /* b(2,6) */
  97. 39520.0/891.0, /* b(3,6) */
  98. -416.0/11.0, /* b(4,6) */
  99. 52.0/27.0, /* b(5,6) */
  100. 0.0, /* b(6,6) */
  101.  
  102. 3015.0/256.0, /* b(1,7) */
  103. -9.0/4.0, /* b(2,7) */
  104. -4219.0/78.0, /* b(3,7) */
  105. 5985.0/128.0, /* b(4,7) */
  106. -539.0/384.0, /* b(5,7) */
  107. 0.0, /* b(6,7) */
  108. 693.0/3328.0, /* b(7,7) */
  109. };
  110.  
  111. #if ERRTERM
  112. double er[INUM] = {
  113. 33.0/640.0, 0.0, -132.0/325.0, 891.0/2240.0, -33.0/320.0,
  114. -73.0/700.0, 891.0/8320.0, 2.0/35.0
  115. };
  116. #endif
  117.  
  118. double a[INUM] = {
  119. 57.0/640.0, 0.0, -16.0/65.0, 1377.0/2240.0, 121.0/320.0, 0.0,
  120. 891.0/8320.0, 2.0/35.0
  121. };
  122. #endif
  123.  
  124.  
  125. #if !(ORDER-7)
  126. /*
  127. Seventh order Runge-Kutta method specified by
  128. Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 7)
  129. */
  130. double al[12] = {
  131. 1.0/4.0, 1.0/12.0, 1.0/8.0, 2.0/5.0, 1.0/2.0, 6.0/7.0, 1.0/7.0,
  132. 2.0/3.0, 2.0/7.0, 1.0, 1.0/3.0, 1.0
  133. };
  134.  
  135. double b[] = {
  136. 1.0/4.0, /* b(1,1) */
  137.  
  138. 5.0/72.0, /* b(1,2) */
  139. 1.0/72.0, /* b(2,2) */
  140.  
  141. 1.0/32.0, /* b(1,3) */
  142. 0.0, /* b(2,3) */
  143. 3.0/32.0, /* b(3,3) */
  144.  
  145. 106.0/125.0, /* b(1,4) */
  146. 0.0, /* b(2,4) */
  147. -408.0/125.0, /* b(3,4) */
  148. 352.0/125.0, /* b(4,4) */
  149.  
  150. 1.0/48.0, /* b(1,5) */
  151. 0.0, /* b(2,5) */
  152. 0.0, /* b(3,5) */
  153. 8.0/33.0, /* b(4,5) */
  154. 125.0/528.0, /* b(5,5) */
  155.  
  156. -1263.0/2401.0, /* b(1,6) */
  157. 0.0, /* b(2,6) */
  158. 0.0, /* b(3,6) */
  159. 39936.0/26411.0, /* b(4,6) */
  160. -64125.0/26411.0, /* b(5,6) */
  161. 5520.0/2401.0, /* b(6,6) */
  162.  
  163. 37.0/392.0, /* b(1,7) */
  164. 0.0, /* b(2,7) */
  165. 0.0, /* b(3,7) */
  166. 0.0, /* b(4,7) */
  167. 1625.0/9408.0, /* b(5,7) */
  168. -2.0/15.0, /* b(6,7) */
  169. 61.0/6720.0, /* b(7,7) */
  170.  
  171. 17176.0/25515.0, /* b(1,8) */
  172. 0.0, /* b(2,8) */
  173. 0.0, /* b(3,8) */
  174. -47104.0/25515.0, /* b(4,8) */
  175. 1325.0/504.0, /* b(5,8) */
  176. -41792.0/25515.0, /* b(6,8) */
  177. 20237.0/145800.0, /* b(7,8) */
  178. 4312.0/6075.0, /* b(8,8) */
  179.  
  180. -23834.0/180075.0, /* b(1,9) */
  181. 0.0, /* b(2,9) */
  182. 0.0, /* b(3,9) */
  183. -77824.0/1980825.0, /* b(4,9) */
  184. -636635.0/633864.0, /* b(5,9) */
  185. 254048.0/300125.0, /* b(6,9) */
  186. -183.0/7000.0, /* b(7,9) */
  187. 8.0/11.0, /* b(8,9) */
  188. -324.0/3773.0, /* b(9,9) */
  189.  
  190. 12733.0/7600.0, /* b(1,10) */
  191. 0.0, /* b(2,10) */
  192. 0.0, /* b(3,10) */
  193. -20032.0/5225.0, /* b(4,10) */
  194. 456485.0/80256.0, /* b(5,10) */
  195. -42599.0/7125.0, /* b(6,10) */
  196. 339227.0/912000.0, /* b(7,10) */
  197. -1029.0/4180.0, /* b(8,10) */
  198. 1701.0/1408.0, /* b(9,10) */
  199. 5145.0/2432.0, /* b(10,10) */
  200.  
  201. -27061.0/204120.0, /* b(1,11) */
  202. 0.0, /* b(2,11) */
  203. 0.0, /* b(3,11) */
  204. 40448.0/280665.0, /* b(4,11) */
  205. -1353775.0/1197504.0, /* b(5,11) */
  206. 17662.0/25515.0, /* b(6,11) */
  207. -71687.0/1166400.0, /* b(7,11) */
  208. 98.0/225.0, /* b(8,11) */
  209. 1.0/16.0, /* b(9,11) */
  210. 3773.0/11664.0, /* b(10,11) */
  211. 0.0, /* 11, 11 */
  212.  
  213. 11203.0/8680.0, /* b(1,12) */
  214. 0.0, /* b(2,12) */
  215. 0.0, /* b(3,12) */
  216. -38144.0/11935.0, /* b(4,12) */
  217. 2354425.0/458304.0, /* b(5,12) */
  218. -84046.0/16275.0, /* b(6,12) */
  219. 673309.0/1636800.0, /* b(7,12) */
  220. 4704.0/8525.0, /* b(8,12) */
  221. 9477.0/10912.0, /* b(9,12) */
  222. -1029.0/992.0, /* b(10,12) */
  223. 0.0, /* b(11,12) */
  224. 729.0/341.0 /* b(12,12) */
  225. };
  226. #if ERRTERM
  227. double er[13] = {
  228. -1.0/480.0, 0.0, 0.0, 0.0, 0.0, -16.0/375.0, -2401.0/528000.0,
  229. 2401.0/132000.0, 243.0/14080.0, -2401.0/19200.0, -19.0/450.0,
  230. 243.0/1760.0, 31.0/720.0
  231. };
  232. #endif
  233. double a[13] = {
  234. 31.0/720.0, 0.0, 0.0, 0.0, 0.0, 16.0/75.0, 16807.0/79200.0,
  235. 16807.0/79200.0, 243.0/1760.0, 0.0, 0.0, 243.0/1760.0, 31.0/720.0
  236. };
  237. #define INUM 13
  238. #define HMAXL 2.0
  239. #define HLIM1 5.0
  240. #endif
  241.  
  242. int inum = INUM;
  243.  
  244. static double *w;
  245. static double *au;
  246. static double *ays;
  247.  
  248. /* Initialize pointers in work array.
  249.  * neq >= neqn in rungek() below.  If neq > neqn, unused space
  250.  * is left for extra variables that are not actually integrated.
  251.  */
  252. rkstart( neq, work )
  253. int neq;
  254. double work[];
  255. {
  256. int m;
  257.  
  258. w = &work[0];
  259. m = neq*inum;
  260. au = &work[m];
  261. m += neq;
  262. ays = &work[m];
  263. }
  264.  
  265.  
  266. /* Make a step.
  267.  */
  268. rungek(  neqn, x, yold, h, ynew, delta )
  269. int neqn; /* number of first order equations */
  270. double x; /* independent variable */
  271. double yold[]; /* initial solution vector at x, size neqn */
  272. double h; /* step to take in independent variable */
  273. double ynew[]; /* output new solution vector at x+h, size neqn */
  274. double delta[]; /* the step, size neqn */
  275. /*double err[];*/ /* estimated error, size neqn */
  276. {
  277. double *pa, *pal, *pb, *per, *pw, *pb0, *pu, *pys;
  278. double xs, usum, esum;
  279. int i, j, m;
  280.  
  281. /*
  282.  * w[]: neqn blocks of size inum
  283.  *      block of size neqn for temp solution vector ys[]
  284.  *      block of size neqn for output function vector u[] of func()
  285.  */
  286. xs = x;
  287.  
  288. pys = ays;
  289. pb = yold;
  290. for( i=0; i<neqn; i++ )
  291.     *pys++ = *pb++;  /* ys[i] = yold[i] */
  292.  
  293. func( xs, ays, au );
  294.  
  295. pw = w;
  296. pu = au;
  297. for( i=0; i<neqn; i++ )
  298.     {
  299.     *pw = *pu++;  /* w[i*inum] = u[i] */
  300.     pw += inum;
  301.     }
  302.  
  303. pb0 = &b[0];
  304. pal = &al[0];
  305. for( j=1; j<inum; j++ )
  306.     {
  307.     xs = x + h * *pal++;
  308.     pys = ays;
  309.     for( i=0; i<neqn; i++ )
  310.         {
  311.         usum = 0.0;
  312.         pw = w + i*inum;
  313.         pb = pb0;
  314.         for( m=0; m<j; m++ )
  315.             usum += (*pb++) * (*pw++);
  316.         *pys++ = yold[i] + h * usum;
  317.         }
  318.     func( xs, ays, au );
  319.     pw = w + j;
  320.     pu = au;
  321.     for( i=0; i<neqn; i++ )
  322.         {
  323.         *pw = *pu++;  /* w[j+i*inum] = u[i] */
  324.         pw += inum;
  325.         }
  326.     pb0 += j;
  327.     }
  328.  
  329. /* evaluate YNEW[i] and the error estimates EST[i] */
  330. pw = w;
  331. for( i=0; i<neqn; i++ )
  332.     {
  333.     pa = &a[0];
  334. #if ERRTERM
  335.     per = &er[0];
  336. #endif
  337.     usum = 0.0;
  338.     esum = 0.0;
  339.     for( j=0; j<inum; j++ )
  340.         {
  341.         if( *pa != 0 )
  342.             {
  343. #if ERRTERM
  344.             esum += (*per++) * (*pw);
  345. #endif
  346.             usum += (*pa++) * (*pw++);
  347.             }
  348.         else
  349.             { /* skip if coefficient = 0 */
  350.             ++pa;
  351. #if ERRTERM
  352.             ++per;
  353. #endif
  354.             ++pw;
  355.             }
  356.         }
  357.     xs = h * usum; /* the step in this coordinate */
  358.     delta[i] = xs;
  359.     ynew[i] = yold[i] + xs;
  360.  
  361. #if ERRTERM
  362.     err[i] = h * esum;
  363. #endif
  364.     }
  365. }
  366.