derkfa dopn43

dopn43() def dopn43 ( info  , n  , f2  , t  , y  , yp  , tout  , tend  , wsave  , iwsave  , mode  = 2, rtol  = 1.0e-10, atol  = 1.0e-10  ) Initial value problem of ordinary differential equations (4(3)-th order Runge-Kutta-Nystrom method) (for second order differential equations) Purpose This program integrates a system of second order ordinary differential equations of the form d2y/dt2 = f(t, y), y = y0, y' = y0' at t = t0 where y and y' are the n vectors and the system consists of n equations. t0, y0 and y'0 are the given initial values of t, y and y', respectively. This program is based on the Runge-Kutta-Nystroem method of order 4(3). See for details of coefficients (RKN4(3)4FM) and interpolant for dense output in the references (1) and (2), respectively. 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. 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] f2 The user supplied subroutine, which calculates the second derivatives of the differential equations, defined as follows. def f2(n, t, y, ypp): Compute d2y/dt2 and store in ypp[i] (for i = 0 to n - 1). where n is the number of equations, and d2y/dt2 (= f2(t, y)) is the computed second 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,out] yp Numpy ndarray (1-dimensional, float, n) [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 8*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 30) 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. (rtol >= 0) (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 and rtol*abs(yp[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. (atol >= 0) (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. y1'' = -y1/R y2'' = -y2/R where R = (y1^2 + y2^2)^(3/2)/(π/4)^2 (y1 = 0.75, y2 = 0, y1' = 0, y2' = (π/4)√(1.25/0.75) at t = 0) def f(n, t, y, ypp): alfa = pi/4 r = (y[0]**2 + y[1]**2)**(3/2)/alfa**2 ypp[0] = -y[0]/r ypp[1] = -y[1]/r def TestDopn43(): n = 2 wsave = np.empty(8*n + 40) iwsave = np.empty(30, dtype=int) alfa = pi/4 ecc = 0.25 t = 0.0 y = np.array([1 - ecc, 0]) yp = np.array([0, alfa*sqrt((1 + ecc)/(1 - ecc))]) tend = 12.0 info = 0 while True: tout = t + 1 t, info = dopn43(info, n, f, t, y, yp, tout, tend, wsave, iwsave) print(t, y, yp) if t >= tend or info < 0 or info > 10: break if info != 0: print('Error: info =', info) Example Results >>> TestDopn43() 1.0 [0.29441755 0.81217852] [-0.7625959 0.47923261] 2.0 [-0.49029978 0.939875 ] [-0.71918028 -0.17238216] 3.0 [-1.05403152 0.5757061 ] [-0.38882952 -0.50909957] 4.0 [-1.25000003e+00 2.61948925e-08] [-3.95363884e-08 -6.08366802e-01] 5.0 [-1.0540316 -0.57570606] [ 0.38882943 -0.5090996 ] 6.0 [-0.49029994 -0.93987504] [ 0.7191802 -0.17238227] 7.0 [ 0.29441735 -0.81217875] [0.76259596 0.47923233] 8.0 [ 7.50000100e-01 -4.82957596e-07] [5.23304832e-07 1.01394459e+00] 9.0 [0.2944181 0.81217828] [-0.76259567 0.47923308] 10.0 [-0.49029923 0.93987522] [-0.71918045 -0.17238178] 11.0 [-1.05403122 0.5757066 ] [-0.38882984 -0.50909938] 12.0 [-1.25000009e+00 6.15458176e-07] [-4.20956957e-07 -6.08366803e-01]