derkfa dopn43

derkfa() def derkfa ( info  , n  , f  , t  , y  , tout  , tend  , wsave  , iwsave  , mode  = 2, rtol  = 1.0e-10, atol  = 1.0e-10  ) Initial value problem of a system of first order ordinary differential equations (5(4)-th order Runge-Kutta-Fehlberg method) Purpose This routine integrates a system of first order ordinary differential equations of the form dy/dt = f(t, y), y = y0 at t = t0 where y is the n vector and the system consists of n equations. t0 and y0 are the given initial values of t and y, respectively. This program is the rewrite of the 5(4)-order Runge-Kutta-Fehrberg program DDERKF (Reference (1)) by adding a dense output function using the interpolation formula of Enright et al. (Reference (2)). Returns (t1, info1) t1 (float): t = tend if the integration has been completed. t = tout or t = (t of the latest step) depending on mode if it is an intermediate return. info1 (int): Return code. = 0: Successful exit. Integration to tend completed. < 0: The (-info)-th argument is invalid. = 1: Returned to provide intermediate result in mode = 1. It is possible to call this program again to forward to the next step. = 2: Returned at t = tout to provide intermediate result in mode = 2. It is possible to call this program again with new tout. = 3: Returned at t = tout to provide intermediate result in mode = 3. It is possible to call this program again with new tout. = 11: (Error) Maximum number of steps exceeded. = 12: (Error) Step size becomes too small. = 13: (Error) Problem is probably stiff. = 14: (Error) Error test failed. Tolerances may be too stringent. Set t = t1 and info = info1 when calling this program again in modes 1, 2 or 3. Parameters [in] info Control code. = 0: Set info = 0 on the initial call to start new problem. All variables will be initialized and the computation will begin. = 1, 2, 3: When returned with info1 = 1, 2 or 3, it is possible to call this program again with setting info = info1 and tout = new tout to continue the integration. [in] n Number of differential equations. (n >= 1) [in] f The user supplied subroutine, which calculates the derivatives of the differential equations, defined as follows. def f(n, t, y, yp): Compute dy/dt and store in yp[i] (for i = 0 to n - 1). where n is the number of equations, and dy/dt (= f(t, y)) is the computed derivative at given t and y. [in] t The current value of the independent variable t (in the first call, equals to the initial value of t = t0). This program integrates from t to tend. Depending on the setting of mode, this program may return at tout or every successful step to provide an intermediate result at t1. [in,out] y Numpy ndarray (1-dimensional, float, n) Dependent variable y. [in] The current value of y at t (In the first call, y at t0 = y0). [out] The value of y (computed numerical solution) at t1. [in] tout Set tout to the point (t) at which an intermediate result is desired in mode = 2 or 3. tout is not referenced in mode = 0 or 1. In mode 2 or 3, y at tout is computed and returned with info1 = 2 or 3 respectively. To continue the integration to get result at new tout, it is possible to call this program again with setting info = info1 and tout = new tout and without changing other variables. t < tout <= tend is required. However, if the integration is in backward direction, tend <= tout < t is required. If tend is reached before tout, the integration is terminated immediately and t1 = tend and info1 = 0 is returned. [out] wsave Numpy ndarray (1-dimensional, float, length 9*n + 40) Data save area (Do not alter while the calculation for the problem starting with info = 0). [out] iwsave Numpy ndarray (1-dimensional, int, length 40) Integer data save area (Do not alter while the calculation for the problem starting with info = 0). [in] mode (Optional) Mode of operation. (default = 2) This program integrates from t to tend. Four modes of operation are provided depending on the timing of returning imtermediate results. When tend is reached, the integration is terminated and info1 = 0 is returned in any of four modes. = 0: Not return until tend. tout is not referenced. = 1: Returns at every successful step (info1 = 1 is returned). tout is not referenced. The end point of the last step is returned in t1. In the final step, the step size is adjusted to fit into tend. = 2: Returns at tout during the ingegration to provide an intermediate result (t1 = tout and info1 = 2 are returned). To continue the integration to get result at new tout, it is possible to call this program again with new tout. In the final step, the step size is adjusted to fit into tout. = 3: Returns at tout during the ingegration to provide an intermediate result (t1 = tout and info1 = 3 are returned). To continue the integration to get result at new tout, it is possible to call this program again with new tout. Different from mode = 2, the value of y at tout is computed by interpolation, and the step size is not adjusted even in the last step before tout. The steps through integration are same with those of mode = 1. [in] rtol (Optional) The relative error tolerance(s) to tell the code how accurately you want the solution to be computed. (default = 1.0e-10) rtol is used together with atol in a local error test at each step. The integration is performed using step sizes which are automatically selected so as to achieve the following equation. abs(local error) <= rtol*abs(y[i]) + atol. Both rtol and atol should not be 0 at the same time. [in] atol (Optional) The absolute error tolerance(s) to tell the code how accurately you want the solution to be computed. (default = 1.0e-10) atol is used together with rtol in a local error test at each step (Refer to rtol above). Reference (1) SLATEC (DEPAC) (2) W. H. Enright, et al.: "Interpolants for Runge-Kutta Formulas", ACM Transactions on Mathematocal Software, Vol.12, No.3. (1986) Example Program Solve the following initial value problem of ordinary differential equations. dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) def f(n, t, y, yp): yp[0] = -2*y[0] + y[1] - cos(t) yp[1] = 2*y[0] - 3*y[1] + 3*cos(t) - sin(t) def TestDerkfa(): n = 2 wsave = np.empty(9*n + 40) iwsave = np.empty(40, dtype=int) t = 0 y = np.array([1.0, 2.0]) tend = 10.0 info = 0 while True: tout = t + 1 t, info = derkfa(info, n, f, t, y, tout, tend, wsave, iwsave) print(t, y) if t >= tend or info < 0 or info > 10: break if info != 0: print('Error: info =', info) Example Results >>> TestDerkf() 1.0 [0.36787944 0.90818175] 2.0 [ 0.13533528 -0.28081155] 3.0 [ 0.04978707 -0.94020543] 4.0 [ 0.01831564 -0.63532798] 5.0 [0.00673795 0.29040013] 6.0 [0.00247875 0.96264904] 7.0 [0.00091188 0.75481414] 8.0 [ 0.00033546 -0.14516457] 9.0 [ 1.23409832e-04 -9.11006852e-01] 10.0 [ 4.53999600e-05 -8.39026129e-01]