arkode

Arguments

Description

arkode computes the solution of real or complex ordinary different equations defined by

The right hand side can be specfied as the sum of two functions where is typically stiff and needs to be integrated by an implicit method, whereas is the non-stiff part of . When both components are specified arkode uses an implicit/explicit method (see the dedicated section below). When the splitting is not necessary a classical explicit (the default) or implicit method can be used (see the available methods below).

The mass matrix M(t) is by default the identity matrix. It can be changed to any non-singular matrix by giving a constant matrix to the mass option or a user function if the matrix is time-dependent. See the mass matrix section below). If the mass matrix depends on solution y then the problem has to be reformulated as a DAE and can be solved by ida.

As arkode shares with cvode most of its features (syntax, specification of user function) please refer to the Description section of cvode help page for basic usage (events handling, solution output, callback). We will give below an example where the solution is computed first with the default explicit method, then with an implicit/explicit method.

Example

In the following example, we solve the Brusselator ODE

The right hand side is computed by the function %SUN_bruss (already defined in the SUNDIALS module):

function dydt=%SUN_bruss(t, y, a, b, ep)
u=y(1); v=y(2); w=y(3);
dydt = [a-(w+1)*u+v*u*u; w*u-v*u*u; (b-w)/ep-w*u];
end

u0 = 3.9; v0 = 1.1; w0 = 2.8; a  = 1.2; b  = 2.5; ep = 1.0e-1;
[t,y] = arkode(list(%SUN_bruss,a,b,ep), [0 10], [u0; v0; w0]);
plot(t,y)

In the example above the default explicit order 4 Runge-Kutta method is used. However, when a very small value is taken for an implicit scheme can be used, here the default order 4 diagonaly implicit Runge-Kutta method:

u0 = 3.9; v0 = 1.1; w0 = 2.8; a  = 1.2; b  = 2.5; ep = 1.0e-5;
[t,y] = arkode(list(%SUN_bruss,a,b,ep), [0 10], [u0; v0; w0], method="DIRK");
plot(t,y)

Using an implicit/explicit method

In the above example we can split the right hand side in two components like this :

function dydt=%SUN_brussI(t, y, b, ep)
    w=y(3);
    dydt = [0; 0; (b-w)/ep];
end
function dydt=%SUN_brussE(t, y, a)
    u=y(1); v=y(2); w=y(3);
    dydt = [a-(w+1)*u+v*u*u; w*u-v*u*u;-w*u];
end

and solve the differential equation with the default order 4 implicit/explicit method. The solver just needs the information about the splitting, which is given by the stiffRhs option:

u0 = 3.9; v0 = 1.1; w0 = 2.8; a  = 1.2; b  = 2.5; ep = 1.0e-5;
jacI = [0 0 0;0 0 0;0 0 -1/ep];
[t,y] = arkode(list(%SUN_brussE,a), [0 10], [u0; v0; w0], stiffRhs = list(%SUN_brussI,b,ep), ... 
               jacobian = jacI);
plot(t,y)

This small example exhibits the interest of the splitting because here the stiff component is linear and we can specify its constant Jacobian to improve the performance of the inner implicit method. Note that when using an ImEx method with a user-provided Jacobian, the latter is the Jacobian of the stiff component only.

Mass matrix

When using arkode a non singular mass matrix can be specified and to this purpose the call to the solver takes the form arkode(f,tspan,y0,mass=massFun) where massFun can be:

• a function with prototype m=massFun(t), for example:

  function m=massFun(t)
      m = [1 0;0 1+t]
  end
  [t,y] = arkode(%SUN_vdp1, [0 10], [0; 2], mass = massFun);

• a matrix, which allows inform the solver that the mass does not depend on time:

  [t,y] = arkode(%SUN_vdp1, [0 10], [0; 2], mass = [1 0;0 2]);

• a string giving the name of a SUNDIALS DLL entrypoint with C prototype

  int sunMass(realtype t, SUNMatrix M, void* user_data, N_Vector tmp1, 
                                        N_Vector tmp2, N_Vector tmp3)

  code=[
      "#include <nvector/nvector_serial.h>"
      "#include <sunmatrix/sunmatrix_dense.h>"
      "int sunRhs(realtype t, N_Vector Y, N_Vector Yd, void *user_data)"    
      "{"
      "double *y = NV_DATA_S(Y);"
      "double *ydot = NV_DATA_S(Yd);"
      "ydot[0] = y[1];"
      "ydot[1] = (1-y[0]*y[0])*y[1]-y[0];"
      "return 0;"
      "}"
      "int sunMass(realtype t, SUNMatrix M, void* user_data, N_Vector tmp1,"
      "            N_Vector tmp2, N_Vector tmp3)"
      "{"
      "double *mass = SM_DATA_D(M);"
      "mass[0] = 1.0; mass[1] = 0.0;"
      "mass[2] = 0.0;  mass[3]= 1.0+t;"
      "return 0;"
      "}"];
  mputl(code,"code.c")
  SUN_Clink(["sunRhs","sunMass"],"code.c",load=%t);
  [t1,y1]=arkode("sunRhs",[0 10],[2;1],mass="sunMass");

Available arkode methods

When there is no additive splitting of the right hand side the default method is the order 4 explicit Zonneveld method (identifiers ERK, ERK_4 or ZONNEVELD_5_3_4). An implicit method can be chosen by using the method option. Note that adding a user Jacobian automatically selects the default implicit method.

If splitting of the right hand side is necessary and the stiff component is given with the stiffRhs option, the default method is a combined order 4 implicit/explicit method (ARK436L2SA_DIRK_6_3_4 and ARK436L2SA_ERK_6_3_4) with identifier ARK or ARK_4.

If you need a custom ERK, DIRK or ARK method see the next section.

Using custom Butcher tableaux

The ERKButcherTab and/or DIRKButcherTab options allow to specify custom Butcher tableaux, under the form of a matrix of size (s+2,s+1) or (s+1,s+1) where s is the number of stages of the method :

A is a matrix of size (s,s), c a column vector of length s, b and d are row vectors of length s, the coefficients of the main and the embedded method, of respective orders q and p. When there is no embedded method (above matrix on the right) the overall matrix has size (s+1,s+1) and adaptive stepping is disabled, hence the option fixedStep must be set to the value of the chosen fixed step. Please note that with a fixed step no error control occurs hence the computed solution may be very inaccurate.

The claimed orders and the structure of A are checked before initialization of the solver: only explicit or diagonally implicit methods are accepted (although the order check allows fully implicit methods). Both options can be used to build a custom additive imEx method, in that case the overall method order can be less than individual ERK/DIRK methods and can be checked in the the arkode solution output.

The explicit Euler and implicit Euler method with fixed step (mandatory in that case) can be specified by:

s1 = arkode(%SUN_vdp1,[0 1],[2;1],ERKButcherTab=[0 0;1 1],fixedStep=0.01);
s2 = arkode(%SUN_vdp1,[0 1],[2;1],DIRKButcherTab=[1 1;1 1],fixedStep=0.01);

The embedded explicit Heun-Euler method can be specified like this,

A = [0 0; 1 0];
c = [0; 1];
b = [0.5 0.5];
d = [1 0];
q = 2;
p = 1;
s = arkode(%SUN_vdp1,[0 10],[2;1],ERKButcherTab=[c A;q b;p d]);

and the implicit trapezoidal rule with implicit Euler embedding has the following tableau

A = [1 0; -1 1];
c = [1; 0];
b = [0.5 0.5];
d = [1 0];
q = 2;
p = 1;
s = arkode(%SUN_vdp1,[0 10],[2;1],DIRKButcherTab=[c A;q b;p d]);

The SUNDIALS module also provides functions %SUN_Tsitouras5, %SUN_DormandPrince6 and %SUN_DormandPrince8 yielding classical tableaux (see references below):

s1 = arkode(%SUN_vdp1,[0 10],[2;1],ERKButcherTab = %SUN_Tsitouras5());
s2 = arkode(%SUN_vdp1,[0 10],[2;1],ERKButcherTab = %SUN_DormandPrince6());
s3 = arkode(%SUN_vdp1,[0 10],[2;1],ERKButcherTab = %SUN_DormandPrince8());

See also

• cvode — SUNDIALS ordinary differential equation solver
• ida — SUNDIALS differential-algebraic equation solver
• User functions — Coding user functions used by SUNDIALS solvers
• Options (ODE and DAE solvers) — Changing the default behavior of solver

Bibliography

A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S. Woodward, "SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers," ACM Transactions on Mathematical Software, 31(3), pp. 363-396, 2005. Also available as LLNL technical report UCRL-JP-200037.

Daniel R. Reynolds and David J. Gardner and Carol S. Woodward and Cody J. Balos, "User Documentation for arkode", 2021, v5.1.1, available online at https://sundials.readthedocs.io/en/latest/arkode.

P.J. Prince and J. R. Dormand, High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics, 7 (1981), pp. 67-75.

Ch. Tsitouras, Runge–Kutta pairs of order 5(4) satisfying only the first column simplifying assumption, Computers & Mathematics with Applications, Volume 62, Issue 2, 2011, Pages 770-775, ISSN 0898-1221, https://doi.org/10.1016/j.camwa.2011.06.002.