derkf derkf_int

derkf() def derkf ( info  , n  , f  , t  , y  , tout  , wsave  , iwsave  , rtol  = 1.0e-10, atol  = 1.0e-10, mode  = 0  ) 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 t0 and y0 are the given initial values of t and y, respectively. y may be a vector if the above is a system of differential equations. This is the code in the package of differential equation solvers DEPAC, consisting of the codes derkf, deabm, and debdf. derkf is a 5(4)-th order Runge-Kutta-Fehlberg code. It is the simplest of the three choices, both algorithmically and in the use of the code. It is primarily designed to solve non-stiff and mildly stiff differential equations when derivative evaluations are not expensive. It should generally not be used to get high accuracy results nor answers at a great many specific points. Because derkf has very low overhead costs, it will usually result in the least expensive integration when solving problems requiring a modest amount of accuracy and having equations that are not costly to evaluate. derkf is a driver for a modification of the code RKF45 written by H. A. Watts and L. F. Shampine. The dense output based on the reference paper below is supported (see mode and example program (2)). Returns (t1, info1) t1 (float): Last value of the independent variable t of the final step (normally equals to tout in interval mode). The solution was successfully advanced to this point. info1 (int): Return code. By examining this code, user can call this routine again with setting t = t1 and info = info1 when info = 1 to 7 as a next action if necessary. = -1: The argument info had an illegal value (info < 0 or info > 7) = -2: The argument n had an illegal value (n < 1) = -3: The argument f is invalid = -4: The argument t had an illegal value (t = tout or t != previous tout on continuation call) = -5: The argument y is invalid = -6: The argument tout had an illegal value (direction of integration is changed) = -7: The argument wsave is invalid = -8: The argument iwsave is invalid = -9: The argument rtol had an illegal value (rtol < 0) = -10: The argument atol had an illegal value (atol < 0) = 1: Successful exit (t = tout). User can reenter with a new tout, which must be different from t. = 2: Interruption in intermediate-output mode (still not reached tout). User can reenter to resume for another step in the direction of tout. = 3: Maximum number of steps (10000) exceeded. User can reenter to resume. Additional 10000 steps will be allowed. = 4: Error tolerances are too stringent. User can reenter with changed tolerances. = 5: Pure relative error test (atol = 0) is impossible because computed solution is zero. User can reenter with positive atol value to continue computation. = 6: Problem is probably stiff. User can reenter to continue, but the other programs such as debdf in DEPAC which handle stiff problems efficiently are recommended. = 7: Natural step size is being restricted by too frequent output. User can reenter to continue, but other programs with dense output feature are recommended. = 8: Infinite loop has been detected. Parameters [in] info Control code. = 0: Set info = 0 on the initial call for the problem (to start new problem). The routine will be initialized and the computation for the new problem will be started. = 1 to 7: When returned from the routine with info1 = 1 to 7, user can reenter with setting info = info1 to continue computation. [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. _CODE def f(n, t, y, yp): yp[i] = computed derivative at t and y (for i = 0 to n-1) _ENDCODE where n is the number of equations, and yp (1-dimensional array with length n) is the derivatives at given t and y (1-dimensional array with length n), i.e. yp[i] = dyi/dt = fi(t, y[0], ..., y[n-1]) (i = 0 to n-1). [in] t Initial value of the independent variable t. This routine integrates from t to tout. The initial point of the integration is to be given, and the last point of the final step will be returned in t1. It is possible to continue the integration with setting t = t1 to get result at new tout. This is the interval mode of operation. t should be equal to the previous tout on continuation call. It is also possible for the routine to return with the solution at each intermediate step on the way to tout. This is the intermediate-output mode of operation. This mode is a good way to proceed if you want to see the behavior of the solution. The mode of operation is specified by the parameter mode. [in,out] y Numpy ndarray (1-dimensional, float, n) [in] Initial values of the dependent variables y at initial t. [out] Computed solution approximation at last t (= tout in interval mode). [in] tout Set tout to the point at which a solution is desired. You can take tout = t, in which case the code will evaluate the derivative of the solution at t and return. Integration either forward in t (tout > t) or backward in t (tout < t) is permitted. It is, however, not allowed to change the direction of integration without restarting. The code advances the solution from t to tout using step sizes which are automatically selected so as to achieve the desired accuracy. If you wish, the code will return with the solution and its derivative following each intermediate step (intermediate-output mode) so that you can monitor them, but you still must provide tout in accord with the basic aim of the code. [out] wsave Numpy ndarray (1-dimensional, float, length 7*n + 20 (9*n + 20 when mode = 2)) Data save area (Do not alter while the calculation for the problem starting with info = 0). [out] iwsave Numpy ndarray (1-dimensional, int, length 20) Integer data save area (Do not alter while the calculation for the problem starting with info = 0). [in,out] rtol (Optional) The relative error tolerance(s) to tell the code how accurately you want the solution to be computed. (default = 1.0e-10) The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error) <= rtol*abs(y[i]) + atol for each component of y (i = 0 to n-1). Setting rtol = 0 yields a pure absolute error test. If you want relative accuracies smaller than about 1.0e-8, you should not ordinarily use derkf. The code deabm in DEPAC obtains stringent accuracies more efficiently. [in,out] atol (Optional) The absolute error tolerance(s) to tell the code how accurately you want the solution to be computed. (default = 1.0e-10) The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error) <= rtol*abs(y[i]) + atol for each component of y (i = 0 to n-1). Setting atol = 0 results in a pure relative error test. If you want relative accuracies smaller than about 1.0e-8, you should not ordinarily use derkf. The code deabm in DEPAC obtains stringent accuracies more efficiently. [in] mode (Optional) Mode of operation. (default = 0) = 0: Return only at tout (interval mode) = 1: Return at every step (intermediate output mode) = 2: Return at every step (intermediate output mode). Dense output is enabled, i.e. derkf_int routine can be used to compute the interpolated values within latest step interval when returned with info1 = 2 (info1 = 1 for last step). (For other values, mode = 0 is assumed.) Reference SLATEC (DEPAC) W H Enright et al. "Interpolants for Runge-Kutta Formulas" ACM Transactions on Mathematical Software Vol.12, No.3, 1986, pp.193-218 Example Program (1) 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 TestDerkf(): n = 2 wsave = np.empty(7*n + 20) iwsave = np.empty(20, dtype=int) t = 0 y = np.array([1.0, 2.0]) tfinal = 10.0 tprint = 1.0 info = 0 while True: tout = t + tprint t, info = derkf(info, n, f, t, y, tout, wsave, iwsave) print(t, y) if t >= tfinal or info != 1: break if info != 1: 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] Example Program (2) Solve the following initial value problem of ordinary differential equations (dense output). 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 TestDerkf_2(): n = 2 wsave = np.empty(9*n + 20) iwsave = np.empty(20, dtype=int) y_temp = np.empty(n) t = 0 y = np.array([1.0, 2.0]) tfinal = 10.0 tprint = 1.0 tout = tprint info = 0 while True: t, info = derkf(info, n, f, t, y, tfinal, wsave, iwsave, mode = 2) if info != 1 and info != 2: print('Error: info =', info) break while t >= tout: derkf_int(n, tout, y_temp, wsave) print(tout, y_temp) tout = tout + tprint if t >= tfinal: break Example Results >>> TestDerkf_2() 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.23409855e-04 -9.11006852e-01] 10.0 [ 4.53999601e-05 -8.39026129e-01]