Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : RUNGEK.C

 
Output of file : RUNGEK.C contained in archive : CEPHES22.ZIP


/* Runge-Kutta numerical integration
* of a system of ordinary differential equations.
*
* Reference:
* Thomas, Benku, "The Runge-Kutta Methods," BYTE 11, #4, April 1986
* This program is adapted from code supplied with that article.
*/

#include "int.h"
#define ORDER 7
#define ERRTERM 0

#if !(ORDER-4)
/*
Fourth order Runge-Kutta method specified by
Felberg, E., COMPUTING, 6(1970)p61-71
*/
#define INUM 6
#define HMAXL 0.5
#define HLIM1 3.0

double al[INUM-1] = {
1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0
};

double b[] = {
1.0/4.0, /* b(1,1) */

3.0/32.0, /* b(1,2) */
9.0/32.0, /* b(2,2) */

1932.0/2197.0, /* b(1,3) */
-7200.0/2197.0, /* b(2,3) */
7296.0/2197.0, /* b(3,3) */

439.0/216.0, /* b(1,4) */
-8.0, /* b(2,4) */
3680.0/513.0, /* b(3,4) */
-845.0/4104.0, /* b(4,4) */

-8.0/27.0, /* b(1,5) */
2.0, /* b(2,5) */
-3544.0/2565.0, /* b(3,5) */
1859.0/4104.0, /* b(4,5) */
-11.0/40.0, /* b(5,5) */
};

#if ERRTERM
double er[INUM] = {
1.0/360.0, 0.0, -384.0/12825.0, -41743.0/1429560.0, 1.0/50.0, 2.0/55.0
};
#endif

double a[INUM] = {
16.0/135.0, 0.0, 6656.0/12825.0, 28561.0/56430.0, -9.0/50.0, 2.0/55.0
};
#endif


#if !(ORDER-5)
/*
Fifth order Runge-Kutta method specified by
Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 5)
*/
#define INUM 8
#define HMAXL 1.0
#define HLIM1 4.0
double al[INUM-1] = {
1.0/18.0, 1.0/6.0, 2.0/9.0, 2.0/3.0, 1.0, 8.0/9.0, 1.0
};

double b[] = {
1.0/18.0, /* b(1,1) */

-1.0/12.0, /* b(1,2) */
1.0/4.0, /* b(2,2) */

-2.0/81.0, /* b(1,3) */
4.0/27.0, /* b(2,3) */
8.0/81.0, /* b(3,3) */

40.0/33.0, /* b(1,4) */
-4.0/11.0, /* b(2,4) */
-56.0/11.0, /* b(3,4) */
54.0/11.0, /* b(4,4) */

-369.0/73.0, /* b(1,5) */
72.0/73.0, /* b(2,5) */
5380.0/219.0, /* b(3,5) */
-12285.0/584.0, /* b(4,5) */
2695.0/1752.0, /* b(5,5) */

-8716.0/891.0, /* b(1,6) */
656.0/297.0, /* b(2,6) */
39520.0/891.0, /* b(3,6) */
-416.0/11.0, /* b(4,6) */
52.0/27.0, /* b(5,6) */
0.0, /* b(6,6) */

3015.0/256.0, /* b(1,7) */
-9.0/4.0, /* b(2,7) */
-4219.0/78.0, /* b(3,7) */
5985.0/128.0, /* b(4,7) */
-539.0/384.0, /* b(5,7) */
0.0, /* b(6,7) */
693.0/3328.0, /* b(7,7) */
};

#if ERRTERM
double er[INUM] = {
33.0/640.0, 0.0, -132.0/325.0, 891.0/2240.0, -33.0/320.0,
-73.0/700.0, 891.0/8320.0, 2.0/35.0
};
#endif

double a[INUM] = {
57.0/640.0, 0.0, -16.0/65.0, 1377.0/2240.0, 121.0/320.0, 0.0,
891.0/8320.0, 2.0/35.0
};
#endif


#if !(ORDER-7)
/*
Seventh order Runge-Kutta method specified by
Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 7)
*/
double al[12] = {
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,
2.0/3.0, 2.0/7.0, 1.0, 1.0/3.0, 1.0
};

double b[] = {
1.0/4.0, /* b(1,1) */

5.0/72.0, /* b(1,2) */
1.0/72.0, /* b(2,2) */

1.0/32.0, /* b(1,3) */
0.0, /* b(2,3) */
3.0/32.0, /* b(3,3) */

106.0/125.0, /* b(1,4) */
0.0, /* b(2,4) */
-408.0/125.0, /* b(3,4) */
352.0/125.0, /* b(4,4) */

1.0/48.0, /* b(1,5) */
0.0, /* b(2,5) */
0.0, /* b(3,5) */
8.0/33.0, /* b(4,5) */
125.0/528.0, /* b(5,5) */

-1263.0/2401.0, /* b(1,6) */
0.0, /* b(2,6) */
0.0, /* b(3,6) */
39936.0/26411.0, /* b(4,6) */
-64125.0/26411.0, /* b(5,6) */
5520.0/2401.0, /* b(6,6) */

37.0/392.0, /* b(1,7) */
0.0, /* b(2,7) */
0.0, /* b(3,7) */
0.0, /* b(4,7) */
1625.0/9408.0, /* b(5,7) */
-2.0/15.0, /* b(6,7) */
61.0/6720.0, /* b(7,7) */

17176.0/25515.0, /* b(1,8) */
0.0, /* b(2,8) */
0.0, /* b(3,8) */
-47104.0/25515.0, /* b(4,8) */
1325.0/504.0, /* b(5,8) */
-41792.0/25515.0, /* b(6,8) */
20237.0/145800.0, /* b(7,8) */
4312.0/6075.0, /* b(8,8) */

-23834.0/180075.0, /* b(1,9) */
0.0, /* b(2,9) */
0.0, /* b(3,9) */
-77824.0/1980825.0, /* b(4,9) */
-636635.0/633864.0, /* b(5,9) */
254048.0/300125.0, /* b(6,9) */
-183.0/7000.0, /* b(7,9) */
8.0/11.0, /* b(8,9) */
-324.0/3773.0, /* b(9,9) */

12733.0/7600.0, /* b(1,10) */
0.0, /* b(2,10) */
0.0, /* b(3,10) */
-20032.0/5225.0, /* b(4,10) */
456485.0/80256.0, /* b(5,10) */
-42599.0/7125.0, /* b(6,10) */
339227.0/912000.0, /* b(7,10) */
-1029.0/4180.0, /* b(8,10) */
1701.0/1408.0, /* b(9,10) */
5145.0/2432.0, /* b(10,10) */

-27061.0/204120.0, /* b(1,11) */
0.0, /* b(2,11) */
0.0, /* b(3,11) */
40448.0/280665.0, /* b(4,11) */
-1353775.0/1197504.0, /* b(5,11) */
17662.0/25515.0, /* b(6,11) */
-71687.0/1166400.0, /* b(7,11) */
98.0/225.0, /* b(8,11) */
1.0/16.0, /* b(9,11) */
3773.0/11664.0, /* b(10,11) */
0.0, /* 11, 11 */

11203.0/8680.0, /* b(1,12) */
0.0, /* b(2,12) */
0.0, /* b(3,12) */
-38144.0/11935.0, /* b(4,12) */
2354425.0/458304.0, /* b(5,12) */
-84046.0/16275.0, /* b(6,12) */
673309.0/1636800.0, /* b(7,12) */
4704.0/8525.0, /* b(8,12) */
9477.0/10912.0, /* b(9,12) */
-1029.0/992.0, /* b(10,12) */
0.0, /* b(11,12) */
729.0/341.0 /* b(12,12) */
};
#if ERRTERM
double er[13] = {
-1.0/480.0, 0.0, 0.0, 0.0, 0.0, -16.0/375.0, -2401.0/528000.0,
2401.0/132000.0, 243.0/14080.0, -2401.0/19200.0, -19.0/450.0,
243.0/1760.0, 31.0/720.0
};
#endif
double a[13] = {
31.0/720.0, 0.0, 0.0, 0.0, 0.0, 16.0/75.0, 16807.0/79200.0,
16807.0/79200.0, 243.0/1760.0, 0.0, 0.0, 243.0/1760.0, 31.0/720.0
};
#define INUM 13
#define HMAXL 2.0
#define HLIM1 5.0
#endif

int inum = INUM;

static double *w;
static double *au;
static double *ays;

/* Initialize pointers in work array.
* neq >= neqn in rungek() below. If neq > neqn, unused space
* is left for extra variables that are not actually integrated.
*/
rkstart( neq, work )
int neq;
double work[];
{
int m;

w = &work[0];
m = neq*inum;
au = &work[m];
m += neq;
ays = &work[m];
}


/* Make a step.
*/
rungek( neqn, x, yold, h, ynew, delta )
int neqn; /* number of first order equations */
double x; /* independent variable */
double yold[]; /* initial solution vector at x, size neqn */
double h; /* step to take in independent variable */
double ynew[]; /* output new solution vector at x+h, size neqn */
double delta[]; /* the step, size neqn */
/*double err[];*/ /* estimated error, size neqn */
{
double *pa, *pal, *pb, *per, *pw, *pb0, *pu, *pys;
double xs, usum, esum;
int i, j, m;

/*
* w[]: neqn blocks of size inum
* block of size neqn for temp solution vector ys[]
* block of size neqn for output function vector u[] of func()
*/
xs = x;

pys = ays;
pb = yold;
for( i=0; i *pys++ = *pb++; /* ys[i] = yold[i] */

func( xs, ays, au );

pw = w;
pu = au;
for( i=0; i {
*pw = *pu++; /* w[i*inum] = u[i] */
pw += inum;
}

pb0 = &b[0];
pal = &al[0];
for( j=1; j {
xs = x + h * *pal++;
pys = ays;
for( i=0; i {
usum = 0.0;
pw = w + i*inum;
pb = pb0;
for( m=0; m usum += (*pb++) * (*pw++);
*pys++ = yold[i] + h * usum;
}
func( xs, ays, au );
pw = w + j;
pu = au;
for( i=0; i {
*pw = *pu++; /* w[j+i*inum] = u[i] */
pw += inum;
}
pb0 += j;
}

/* evaluate YNEW[i] and the error estimates EST[i] */
pw = w;
for( i=0; i {
pa = &a[0];
#if ERRTERM
per = &er[0];
#endif
usum = 0.0;
esum = 0.0;
for( j=0; j {
if( *pa != 0 )
{
#if ERRTERM
esum += (*per++) * (*pw);
#endif
usum += (*pa++) * (*pw++);
}
else
{ /* skip if coefficient = 0 */
++pa;
#if ERRTERM
++per;
#endif
++pw;
}
}
xs = h * usum; /* the step in this coordinate */
delta[i] = xs;
ynew[i] = yold[i] + xs;

#if ERRTERM
err[i] = h * esum;
#endif
}
}


  3 Responses to “Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : RUNGEK.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/