[Search for users] [Overall Top Noters] [List of all Conferences] [Download this site]

Conference rusure::math

Title:Mathematics at DEC
Moderator:RUSURE::EDP
Created:Mon Feb 03 1986
Last Modified:Fri Jun 06 1997
Last Successful Update:Fri Jun 06 1997
Number of topics:2083
Total number of notes:14613

1256.0. "ODE solver" by HDLITE::ALJAAR () Thu Jun 21 1990 12:03

Hi;

	Can anyone point me to an ordinary differential equation solver that
runs under VMS? Also, if anyone knows of a MAC-based ODE solver, please send
me info.

	Thanks..

Robi.
T.RTitleUserPersonal
Name		DateLines
1256.1Numerical Recipes strikes again!ALLVAX::JROTHIt's a bush recording...Thu Jun 21 1990 22:1615
�	Can anyone point me to an ordinary differential equation solver that
� runs under VMS? Also, if anyone knows of a MAC-based ODE solver, please send
� me info.

    If you have a normal, non-stiff problem, then use a Runge Kutta
    algorithm, such as RK45, out of the Numerical Recipes book.  It's really
    easy to do - the integrator routine is often smaller than the function
    you're integrating! (well, almost.)

    I have some fairly general routines online and could mail them if
    you wish, but it would almost be easier to use the Recipes routine
    than the general ones.  I'd post them but they are probably a bit long
    to fill up the notes file with.

    - Jim
1256.2Thanks...HDLITE::ALJAARFri Jun 22 1990 16:520
1256.3RUNGE-KUTA & ADAMS-MOULTON method example in CSUBWAY::TJIONASGeorge Tjionas - NY PSSMon Jun 25 1990 17:30441
$ type diff_equation.c

#include <stdio.h>
#include <math.h>

double float    xfirst  = 0.00;    /* lower bound of x */
double float    xlast   = 1.00;    /* upper bound of x */
double float    yfirst  = 0.40;    /* initial value of y */
double float    h       = 0.10;    /* increment value */
int             i;                 /* looping index */

double float    x[100];             /* for x[i] */

double float    y[4];               /* for RUNGE-KUTTA y[i] */
double float    f[4];               /* for RUNGE-KUTTA f(x[i], y[i]) */
double float    yp[100];            /* for ADAMS-MOULTON predictor yc[i] */
double float    yc[100];            /* for ADAMS-MOULTON corrector yc[i] */
double float    fp[100];            /* for f-predictor f(xp[i], yp[i]) */
double float    fc[100];            /* for f-corrector f(xc[i], yc[i]) */

double float    k1      = 0.0;     /* k1 = h*f[100] = h*f(x[i], y[i]) */
double float    x1      = 0.0;     /* x1 = x[i] + h/2  */
double float    y1      = 0.0;     /* y1 = y[i] + k1/2 */
double float    f1      = 0.0;     /* f1 = f(x[i]+h/2, y[i]+k1/2) */

double float    k2      = 0.0;     /* k2 = h*f(x[i]+h/2, y[i]+k1/2) */
                                   /* k2 = h*f(x1, y1) */
double float    x2      = 0.0;     /* x2 = x[i] + h/2 = x1  */
double float    y2      = 0.0;     /* y2 = y[i] + k2/2 */
double float    f2      = 0.0;     /* f2 = f(x[i]+h/2, y[i]+k2/2) */
                                   /* f2 = f(x1, y2) */

double float    k3      = 0.0;     /* k3 = h*f(x[i]+h/2, y[i]+k2/2) */
                                   /* k3 = h*f(x1, y2) */
double float    x3      = 0.0;     /* x3 = x[i] + h  */
double float    y3      = 0.0;     /* y3 = y[i] + k3 */
double float    f3      = 0.0;     /* f3 = f(x[i]+h, y[i]+k3) */
                                   /* f3 = f(x3, y3) */

double float    k4      = 0.0;     /* k4 = h*f(x[i]+h, y[i]+k3) */
                                   /* k4 = h*f(x3, y3) */

double float    y4      = 0.0;     /* y4 = y[i]+(k1+2*k2+2*k3+k4)/6 = y[i+1] */
double float    x4      = 0.0;     /* x4 = x[i] + h/2 = x[i+1]  */


main()
{
double float fxy_equation();


/* initialize RUNGE-KUTTA arrays */
for (i=0; i<=3; i++)
    {
	y[i] = 0.0;             /* for y[i] */
	f[i] = 0.0;             /* for f(x[i], y[i]) */
    }; /* end of for loop */

/* initialize ADAMS-MOULTON arrays */
for (i=0; i<=99; i++)
    {
	x[i]  = 0.0;            /* for x[i] */
	yp[i] = 0.0;            /* for y[i] */
	yc[i] = 0.0;            /* for y[i] */
	fp[i] = 0.0;            /* for f(x[i], y[i]) */
	fc[i] = 0.0;            /* for f(x[i], y[i]) */
    }; /* end of for loop */

/* set initial values */
	x[0] = xfirst;
        y[0] = yfirst;


printf("\n  i     h      Xi           Yi\n");        /* header */
printf(   "----  ----    -----    -------------\n");  /* underline */
/* RUNGE-KUTTA method */
for (i=0; i<=3; i++)
    {
	if ( x[i] > xlast )
		break;	/* termination criterion */

	printf(" %2d   %1.2f    %1.3f    %2.10f\n", i, h, x[i], y[i]);

	f[i]  = fxy_equation(x[i], y[i]);
	yp[i] = yc[i] = y[i];                /* save for ADAMS-MOULTON */
	fp[i] = fc[i] = f[i];                /* save for ADAMS-MOULTON */

	k1 = h*f[i];
	x1 = x[i] + h/2;
	y1 = y[i] + k1/2;
	f1 = fxy_equation(x[i]+h/2, y[i]+k1/2);

	k2 = h*fxy_equation(x[i]+h/2, y[i]+k1/2);/* k2=h*fxy_equation(x1,y1) */
	x2 = x[i] + h/2;
	y2 = y[i] + k2/2;
	f2 = fxy_equation(x[i]+h/2, y[i]+k2/2);  /* f2=fxy_equation(x1,y2) */

	k3 = h*fxy_equation(x[i]+h/2, y[i]+k2/2);/* k3=h*fxy_equation(x1,y2) */
	x3 = x[i] + h;
	y3 = y[i] + k3;
	f3 = fxy_equation(x[i]+h, y[i]+k3);       /* f3=fxy_equation(x3,y3) */

	k4 = h*fxy_equation(x[i]+h, y[i]+k3);     /* k4=h*fxy_equation(x3,y3) */

	y4 = y[i]+(k1+2*k2+2*k3+k4)/6;

	y[i+1] = y4;
	x[i+1] = x3;

    }; /* end of for loop */

if ( x[i] > xlast )
	return;	/* termination criterion - x exceeds upper bound */

/* ADAMS-MOULTON method */
/* Note: At this point, the i was left as i = 4 from the for loop */
i = 3;   /* set it back since the calculations in the formula refer to i+1 */
while ( x[i] <= xlast && i <= 99 )    /* while x is within range */
    {
	x[i+1] = x[i] + h;

	/* predictor */
	yp[i+1]=yp[i] + (h/24)*(55*fc[i] - 59*fc[i-1] + 37*fc[i-2] - 9*fc[i-3]);
	fp[i+1]=fxy_equation(x[i+1], yp[i+1]);

	/* corrector */
	yc[i+1]=yc[i] + (h/24)*(9*fp[i+1] + 19*fc[i] - 5*fc[i-1] + fc[i-2]);
	fc[i+1]=fxy_equation(x[i+1], yc[i+1]);

	i += 1;
	printf(" %2d   %1.2f    %1.3f    %2.10f\n", i, h, x[i], yc[i]);
    }; /* end of for loop */

} /* end main() */


double float fxy_equation( double float x, double float y )
{
	/* This is my differential function */
        /*                       2 */
	/*  f(x, y) = (1 + 2*x)*y  */

	/* the math pow() function is the power of base x to exponent y */

        return( (1 + 2*x)*(pow(y, 2.0)) );

} /* end of double float fxy_equation() function */

$ run diff_equation

  i     h      Xi           Yi
----  ----    -----    -------------
  0   0.10    0.000    0.4000000000
  1   0.10    0.100    0.4184100695
  2   0.10    0.200    0.4424779459
  3   0.10    0.300    0.4739337793
  4   0.10    0.400    0.5154729556
  5   0.10    0.500    0.5714393601
  6   0.10    0.600    0.6493330468
  7   0.10    0.700    0.7631971307
  8   0.10    0.800    0.9425913247
  9   0.10    0.900    1.2615925379
 10   0.10    1.000    1.9671322970



$ type diff_equation1.c

#include <stdio.h>
#include <math.h>

double float    h;                 /* increment value */
int             i;                 /* looping index */

double float    x[10];             /* for x[i] */
double float    y[10];             /* for y[i] */
double float    f[10];             /* for f(x[i], y[i]) */

double float    k1      = 0.0;     /* k1 = h*f(x[i], y[i]) */
double float    x1      = 0.0;     /* x1 = x[i] + h/2  */
double float    y1      = 0.0;     /* y1 = y[i] + k1/2 */
double float    f1      = 0.0;     /* f1 = f(x[i]+h/2, y[i]+k1/2) */

double float    k2      = 0.0;     /* k2 = h*f(x[i]+h/2, y[i]+k1/2) */
                                   /* k2 = h*f(x1, y1) */
double float    x2      = 0.0;     /* x2 = x[i] + h/2 = x1  */
double float    y2      = 0.0;     /* y2 = y[i] + k2/2 */
double float    f2      = 0.0;     /* f2 = f(x[i]+h/2, y[i]+k2/2) */
                                   /* f2 = f(x1, y2) */

double float    k3      = 0.0;     /* k3 = h*f(x[i]+h/2, y[i]+k2/2) */
                                   /* k3 = h*f(x1, y2) */
double float    x3      = 0.0;     /* x3 = x[i] + h  */
double float    y3      = 0.0;     /* y3 = y[i] + k3 */
double float    f3      = 0.0;     /* f3 = f(x[i]+h, y[i]+k3) */
                                   /* f3 = f(x3, y3) */

double float    k4      = 0.0;     /* k4 = h*f(x[i]+h, y[i]+k3) */
                                   /* k4 = h*f(x3, y3) */

double float    y4      = 0.0;     /* y4 = y[i]+(k1+2*k2+2*k3+k4)/6 = y[i+1] */
double float    x4      = 0.0;     /* x4 = x[i] + h/2 = x[i+1]  */


/* RUNGE_KUTTA method */
main()
{
double float fxy_equation();


/* initialize arrays */
for (i=0; i<=3; i++)
    {
	x[i] = 0.0;             /* for x[i] */
	y[i] = 0.0;             /* for y[i] */
	f[i] = 0.0;             /* for f(x[i], y[i]) */
    }; /* end of for loop */

/* set initial values */
	h    = 0.10;
	x[0] = 1.0;
        y[0] = 0.0;

/* do calculations */
printf("\n  i     h      Xi          Yi\n");        /* header */
printf(   "----  ----    -----    ------------\n");  /* underline */
for (i=0; i<=3; i++)
    {
	printf(" %2d   %1.2f    %1.3f    %1.10f\n", i, h, x[i], y[i]);

	k1 = h*fxy_equation(x[i], y[i]);
	x1 = x[i] + h/2;
	y1 = y[i] + k1/2;
	f1 = fxy_equation(x[i]+h/2, y[i]+k1/2);

	k2 = h*fxy_equation(x[i]+h/2, y[i]+k1/2); /* k2=h*fxy_equation(x1,y1) */
	x2 = x[i] + h/2;
	y2 = y[i] + k2/2;
	f2 = fxy_equation(x[i]+h/2, y[i]+k2/2);   /* f2=fxy_equation(x1,y2) */

	k3 = h*fxy_equation(x[i]+h/2, y[i]+k2/2); /* k3=h*fxy_equation(x1,y2) */
	x3 = x[i] + h;
	y3 = y[i] + k3;
	f3 = fxy_equation(x[i]+h, y[i]+k3);       /* f3=fxy_equation(x3,y3) */

	k4 = h*fxy_equation(x[i]+h, y[i]+k3);     /* k4=h*fxy_equation(x3,y3) */

	y4 = y[i]+(k1+2*k2+2*k3+k4)/6;

	y[i+1] = y4;
	x[i+1] = x3;
    }; /* end of for loop */

} /* end main() */


double float fxy_equation( double float x, double float y )
{
/* This is my differential function */
        /*             y */
	/*  f(x, y) = x  */

        return( pow(x, y) );   /* power of base x to  exponent y */

} /* end of double float fxy_equation() function */

$ run diff_equation1

  i     h      Xi          Yi
----  ----    -----    ------------
  0   0.10    1.000    0.0000000000
  1   0.10    1.100    0.1003230339
  2   0.10    1.200    0.2025326744
  3   0.10    1.300    0.3084666131



$ type diff_equation2.c

#include <stdio.h>
#include <math.h>

double float    xfirst  = 1.0;     /* lower bound of x */
double float    xlast   = 2.0;     /* upper bound of x */
double float    yfirst  = 0.0;     /* initial value of y */
double float    h       = 0.10;    /* increment value */
int             i;                 /* looping index */

double float    x[100];             /* for x[i] */

double float    y[4];               /* for RUNGE-KUTTA y[i] */
double float    f[4];               /* for RUNGE-KUTTA f(x[i], y[i]) */
double float    yp[100];            /* for ADAMS-MOULTON predictor yc[i] */
double float    yc[100];            /* for ADAMS-MOULTON corrector yc[i] */
double float    fp[100];            /* for f-predictor f(xp[i], yp[i]) */
double float    fc[100];            /* for f-corrector f(xc[i], yc[i]) */

double float    k1      = 0.0;     /* k1 = h*f[100] = h*f(x[i], y[i]) */
double float    x1      = 0.0;     /* x1 = x[i] + h/2  */
double float    y1      = 0.0;     /* y1 = y[i] + k1/2 */
double float    f1      = 0.0;     /* f1 = f(x[i]+h/2, y[i]+k1/2) */

double float    k2      = 0.0;     /* k2 = h*f(x[i]+h/2, y[i]+k1/2) */
                                   /* k2 = h*f(x1, y1) */
double float    x2      = 0.0;     /* x2 = x[i] + h/2 = x1  */
double float    y2      = 0.0;     /* y2 = y[i] + k2/2 */
double float    f2      = 0.0;     /* f2 = f(x[i]+h/2, y[i]+k2/2) */
                                   /* f2 = f(x1, y2) */

double float    k3      = 0.0;     /* k3 = h*f(x[i]+h/2, y[i]+k2/2) */
                                   /* k3 = h*f(x1, y2) */
double float    x3      = 0.0;     /* x3 = x[i] + h  */
double float    y3      = 0.0;     /* y3 = y[i] + k3 */
double float    f3      = 0.0;     /* f3 = f(x[i]+h, y[i]+k3) */
                                   /* f3 = f(x3, y3) */

double float    k4      = 0.0;     /* k4 = h*f(x[i]+h, y[i]+k3) */
                                   /* k4 = h*f(x3, y3) */

double float    y4      = 0.0;     /* y4 = y[i]+(k1+2*k2+2*k3+k4)/6 = y[i+1] */
double float    x4      = 0.0;     /* x4 = x[i] + h/2 = x[i+1]  */


main()
{
double float fxy_equation();


/* initialize RUNGE-KUTTA arrays */
for (i=0; i<=3; i++)
    {
	y[i] = 0.0;             /* for y[i] */
	f[i] = 0.0;             /* for f(x[i], y[i]) */
    }; /* end of for loop */

/* initialize ADAMS-MOULTON arrays */
for (i=0; i<=99; i++)
    {
	x[i]  = 0.0;            /* for x[i] */
	yp[i] = 0.0;            /* for y[i] */
	yc[i] = 0.0;            /* for y[i] */
	fp[i] = 0.0;            /* for f(x[i], y[i]) */
	fc[i] = 0.0;            /* for f(x[i], y[i]) */
    }; /* end of for loop */

/* set initial values */
	x[0] = xfirst;
        y[0] = yfirst;


printf("\n  i     h      Xi           Yi\n");        /* header */
printf(   "----  ----    -----    -------------\n");  /* underline */
/* RUNGE-KUTTA method */
for (i=0; i<=3; i++)
    {
	if ( x[i] > xlast )
		break;	/* termination criterion */

	printf(" %2d   %1.2f    %1.3f    %2.10f\n", i, h, x[i], y[i]);

	f[i]  = fxy_equation(x[i], y[i]);
	yp[i] = yc[i] = y[i];                /* save for ADAMS-MOULTON */
	fp[i] = fc[i] = f[i];                /* save for ADAMS-MOULTON */

	k1 = h*f[i];
	x1 = x[i] + h/2;
	y1 = y[i] + k1/2;
	f1 = fxy_equation(x[i]+h/2, y[i]+k1/2);

	k2 = h*fxy_equation(x[i]+h/2, y[i]+k1/2);/* k2=h*fxy_equation(x1,y1) */
	x2 = x[i] + h/2;
	y2 = y[i] + k2/2;
	f2 = fxy_equation(x[i]+h/2, y[i]+k2/2);  /* f2=fxy_equation(x1,y2) */

	k3 = h*fxy_equation(x[i]+h/2, y[i]+k2/2);/* k3=h*fxy_equation(x1,y2) */
	x3 = x[i] + h;
	y3 = y[i] + k3;
	f3 = fxy_equation(x[i]+h, y[i]+k3);       /* f3=fxy_equation(x3,y3) */

	k4 = h*fxy_equation(x[i]+h, y[i]+k3);     /* k4=h*fxy_equation(x3,y3) */

	y4 = y[i]+(k1+2*k2+2*k3+k4)/6;

	y[i+1] = y4;
	x[i+1] = x3;

    }; /* end of for loop */

if ( x[i] > xlast )
	return;	/* termination criterion - x exceeds upper bound */

/* ADAMS-MOULTON method */
/* Note: At this point, the i was left as i = 4 from the for loop */
i = 3;   /* set it back since the calculations in the formula refer to i+1 */
while ( x[i] <= xlast && i <= 99 )    /* while x is within range */
    {
	x[i+1] = x[i] + h;

	/* predictor */
	yp[i+1]=yp[i] + (h/24)*(55*fc[i] - 59*fc[i-1] + 37*fc[i-2] - 9*fc[i-3]);
	fp[i+1]=fxy_equation(x[i+1], yp[i+1]);

	/* corrector */
	yc[i+1]=yc[i] + (h/24)*(9*fp[i+1] + 19*fc[i] - 5*fc[i-1] + fc[i-2]);
	fc[i+1]=fxy_equation(x[i+1], yc[i+1]);

	i += 1;
	printf(" %2d   %1.2f    %1.3f    %2.10f\n", i, h, x[i], yc[i]);
    }; /* end of for loop */

} /* end main() */


double float fxy_equation( double float x, double float y )
{
double float   e = 2.718281828;
	/* This is my differential function */
        /*                     2  x */
	/*  f(x, y) = 2*y/x + x *e  */

	/* the math pow() function is the power of base x to exponent y */

        return( (2*y)/x + pow(x, 2.0)*pow(e, x) );

} /* end of double float fxy_equation() function */

$ run diff_equation2

  i     h      Xi           Yi
----  ----    -----    -------------
  0   0.10    1.000    0.0000000000
  1   0.10    1.100    0.3459102872
  2   0.10    1.200    0.8666216926
  3   0.10    1.300    1.6071813473
  4   0.10    1.400    2.6203319407
  5   0.10    1.500    3.9676301772
  6   0.10    1.600    5.7208998722
  7   0.10    1.700    7.9637681672
  8   0.10    1.800    10.7934560940
  9   0.10    1.900    14.3228284861
 10   0.10    2.000    18.6827365614
1256.4Thanks!HDLITE::ALJAARTue Jun 26 1990 16:380
1256.5try MapleCIVAGE::LYNNLynn Yarbrough @WNP DTN 427-5663Mon Jul 02 1990 12:445
For some equations, you can get an explicit, as opposed to a strictly 
numerical, solution with something like MAPLE. It might help in getting an 
insight to the form of the general solution.

Lynn Yarbrough