C 代码 求解常微分方程组（ODE）.rar

# include <float.h> # include <math.h> # include <stdio.h> # include <stdlib.h> # include <time.h> # include "ode.h" /******************************************************************************/ void de ( void f ( double t, double y[], double yp[] ), int neqn, double y[], double *t, double tout, double relerr, double abserr, int *iflag, double yy[], double wt[], double p[], double yp[], double ypout[], double phi[], double alpha[], double beta[], double sig[], double v[], double w[], double g[], int *phase1, double psi[], double *x, double *h, double *hold, int *start, double *told, double *delsgn, int *ns, int *nornd, int *k, int *kold, int *isnold ) /******************************************************************************/ /* Purpose: DE carries out the ODE solution algorithm. Discussion: ODE merely allocates storage for DE, to relieve the user of the inconvenience of a long call list. Consequently, DE is used as described in the comments for ODE. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 February 2012 Author: Original FORTRAN77 version by Lawrence Shampine, Marilyn Gordon. C version by John Burkardt. Reference: Lawrence Shampine, Marilyn Gordon, Computer Solution of Ordinary Differential Equations: The Initial Value Problem, Freeman, 1975, ISBN: 0716704617, LC: QA372.S416. Parameters: Input, void F ( double T, double Y[], double YP[] ), the user-supplied function which accepts input values T and Y[], evaluates the right hand sides of the ODE, and stores the result in YP[]. Input, int NEQN, the number of equations. Input/output, double Y[NEQN], the current solution. Input/output, double *T, the current value of the independent variable. Input, double TOUT, the desired value of T on output. Input, double RELERR, ABSERR, the relative and absolute error tolerances. At each step, the code requires abs ( local error ) <= abs ( Y ) * RELERR + ABSERR for each component of the local error and solution vectors. Input/output, int *IFLAG, indicates the status of integration. On input, IFLAG is normally 1 (or -1 in the special case where TOUT is not to be exceeded.) On normal output, IFLAG is 2. Other output values are: * 3, integration did not reach TOUT because the error tolerances were too small. But RELERR and ABSERR were increased appropriately for continuing; * 4, integration did not reach TOUT because more than 500 steps were taken; * 5, integration did not reach TOUT because the equations appear to be stiff; * 6, invalid input parameters (fatal error). The value of IFLAG is returned negative when the input value is negative and the integration does not reach TOUT. Workspace, double YY[NEQN], used to hold old solution data. Input, double WT[NEQN], the error weight vector. Workspace, double P[NEQN]. Workspace, double YP[NEQN], used to hold values of the solution derivative. Workspace, double YPOUT[NEQN], used to hold values of the solution derivative. Workspace, double PHI[NEQN*16], contains divided difference information about the polynomial interpolant to the solution. Workspace, double ALPHA[12], BETA[12], SIG[13]. Workspace, double V[12], W[12], G[13]. Input/output, int *PHASE1, indicates whether the program is in the first phase, when it always wants to increase the ODE method order. Workspace, double PSI[12], contains information about the polynomial interpolant to the solution. Input/output, double *X, a "working copy" of T, the current value of the independent variable, which is adjusted as the code attempts to take a step. Input/output, double *H, the current stepsize. Input/output, double *HOLD, the last successful stepsize. Input/output, int *START, is TRUE on input for the first step. The program initializes data, and sets START to FALSE. Input/output, double *TOLD, the previous value of T. Input/output, double *DELSGN, the sign (+1 or -1) of TOUT - T. Input/output, int *NS, the number of steps taken with stepsize H. Input/output, int *NORND, ? Input/output, int *K, the order of the current ODE method. Input/output, int *KOLD, the order of the ODE method on the previous step. Input/output, int *ISNOLD, the previous value of ISN, the sign of IFLAG. Local parameters: Local, integer MAXNUM, the maximum number of steps allowed in one call to DE. */ { double absdel; double abseps; int crash; double del; double eps; double fouru; int isn; int kle4; int l; const int maxnum = 500; int nostep; double releps; int stiff; double tend; /* Test for improper parameters. */ fouru = 4.0 * DBL_EPSILON; if ( neqn < 1 ) { *iflag = 6; fprintf ( stderr, "\n" ); fprintf ( stderr, "DE - Fatal error!\n" ); fprintf ( stderr, " NEQN < 1.\n" ); exit ( 1 ); } if ( *t == tout ) { *iflag = 6; fprintf ( stderr, "\n" ); fprintf ( stderr, "DE - Fatal error!\n" ); fprintf ( stderr, " T = TOUT.\n" ); exit ( 1 ); } if ( relerr < 0.0 || abserr < 0.0 ) { *iflag = 6; fprintf ( stderr, "\n" ); fprintf ( stderr, "DE - Fatal error!\n" ); fprintf ( stderr, " RELERR < 0 or ABSERR < 0.\n" ); exit ( 1 ); } eps = fmax ( relerr, abserr ); if ( eps <= 0.0 ) { *iflag = 6; fprintf ( stderr, "\n" ); fprintf ( stderr, "DE - Fatal error!\n" ); fprintf ( stderr, " max ( RELERR, ABSERR ) <= 0.\n" ); exit ( 1 ); } if ( *iflag == 0 ) { *iflag = 6; fprintf ( stderr, "\n" ); fprintf ( stderr, "DE - Fatal error!\n" ); fprintf ( stderr, " IFLAG = 0 on input.\n" ); exit ( 1 ); } isn = i4_sign ( *iflag ); *iflag = abs ( *iflag ); if ( *iflag != 1 ) { if ( *t != *told ) { *iflag = 6; fprintf ( stderr, "\n" ); fprintf ( stderr, "DE - Fatal error!\n" ); fprintf ( stderr, " IFLAG is not 1, and T is not equal to TOLD.\n" ); exit ( 1 ); } if ( *iflag < 2 || 5 < *iflag ) { *iflag = 6; return; } } /* On each call set interval of integration and counter for number of steps. Adjust input error tolerances to define weight vector for subroutine STEP. */ del = tout - *t; absdel = fabs ( del ); if ( isn < 0 ) { tend = tout; } else { tend = *t + 10.0 * del; } nostep = 0; kle4 = 0; stiff = 0; releps = relerr / eps; abseps = abserr / eps; /* On start and restart, also set work variables X and YY(*), store the direction of integration, and initialize the step size. */ if ( *iflag == 1 || *isnold < 0 || *delsgn * del <= 0.0 ) { *start = 1; *x = *t; for ( l = 1; l <= neqn; l++ ) { yy[l-1] = y[l-1]; } *delsgn = r8_sign ( del ); *h = fmax ( fabs ( tout - *x ), fouru * fabs ( *x ) ) * r8_sign ( tout - *x ); } /* If already past the output point, then interpolate and return. */ for ( ; ; ) { if ( absdel <= fabs ( *x - *t ) ) { intrp ( *x, yy, tout, y, ypout, neqn, *kold, phi, psi ); *iflag = 2; *t = tout; *told = *t; *isnold = isn; break; } /* If we cannot go past the output point, and we are sufficiently close to it, then extrapolate and return. */ if ( isn <= 0 && fabs ( tout - *x ) < fouru * fabs ( *x ) ) { *h = tout - *x; f ( *x, yy, yp ); for ( l = 1; l <= neqn; l++ ) { y[l-1] = yy[l-1] + *h * yp[l-1]; } *iflag = 2; *t = tout; *told = *t; *isnold = isn; break; } /* Test for too many steps. */ if ( maxnum <= nostep ) { *iflag = isn * 4; if ( stiff )