Differential Equations
Mathwrist C++ API Series

Description:

It returns the dimension size $n$ of the ODE system, which is the number of unknown variables $y(t)$ to be solved simultaneously.

Description:

It returns the ODE system start time $t_0$ and end time $T$ respectively.

Description:

It computes the derivative $y^{\prime }(t)$ as of time $t$, given the values of $y(t)$. Both parameters _y and _fty are pointers to $n$-element buffer that are managed by NPL ODE solvers. Pointer _y stores the values of $y(t)$. On return, parameter _fty holds the values of $f(t,y(t))$.

Description:

Return true if the ODE is homogeneous, $s(t)=0$.

Description:

$A(t)$ very often is a sparse matrix. Users may choose a compact storage to represent $A(t)$ and specify the total number of elements of $A(t)$ in their compact form. NPL solvers then use this information to allocate enough memory to store $A(t)$. If $A(t)$ is constant, a derived class may return the size of $A(t)$ to be 0 and cache the constant matrix $A(t)$ internally.

Description:

This function evaluates $A(t)$ and $s(t)$ as of time $t$ and stores the result to provided buffer pointers _A and _s.

Parameters _A points to a memory buffer managed by NPL solvers. _A is allocated to store the number of elements specified by size_A() function. When size_A() returns 0, pointer _A will be a null pointer. In which case, NPL assumes $A(t)$ is handled internally in a derived class.

When the ODE is homogeneous, parameter _s is passed as a null pointer and can be ignored. Otherwise, it points to a memory buffer managed by NPL solvers. It should store $n$-element of source function $s(t)$ values on return.

Description:

This function asks the concrete linear ODE class to compute $y^{\prime }(t)=A(t)y(t)+s(t)$ with given $A(t)$, $y(t)$ and $s(t)$. Parameters _A and _s are same as that of the update(t, _A, _s) function. Parameter _y points to a buffer that stores $n$-element of $y(t)$ values. Parameter _yp is pointer to $n$-element buffer that stores the result $y^{\prime }(t)$ on return.

Description:

This function requires a concrete linear ODE derived class to solve linear system (2) using its compact storage $A(t)$ that is passed from parameter _A. Parameter t is the time variable value. Parameter h is ODE solver’s time marching step length. Parameter _y holds the $n$-element right hand side vector $𝐫$ on call. It should store the solution vector $𝐲$ on return.

Description

This ODE_IDeC::Method enumeration list is a member of class ODE_IDeC It gives users different combination choices of base methods and IDeC methods. A base method starts the process with an first-order or second order approximation of $y(t)$. IDeC methods then raise the accuracy order at each error correction iteration.

Description:

Construct an IDeC solver object given a method of choice, degree of interpolation/approximation polynomial and number of error correction iterations.

Parameter m is one of the ODE_IDeC::Method enumeration item. It specifies what base method and IDeC method are used.

Parameter p specifies the degree of interpolation polynomial used in class IDeC method or the degree of Chebyshev approximation polynomial in spectral IDeC. REQUIRED: in classic IDeC and in spectral IDeC.

Parameter q specifies the number of error correction iterations. $q=0$ implies no error correction is applied, in which case the approximation result from the base method is immediately returned. REQUIRED: .

Description:

Solve an ODE problem by iterative deferred correction method.

Parameter o is an object of concrete ODE problems supplied by the user. It is derived from the ODE base class for problems to be solved with explicit methods, or derived from ODE_Linear class if implicit methods are desired.

Parameter N is the suggested number of time discretization steps. For problems to be solved with explicit methods, time marching step size needs be fine enough to retain numerical stability.

Parameter y0 is a $n$-elment buffer managed by the caller to provide initial values of $y(0)$.

Parameter yT is a $n$-element buffer managed by the caller. It stores the numerical solution $y(T)$ on return.

Description:

The nested class PDEParabolic1D::BoundaryCondition represents boundary conditions in the general Robin boundary condition form. For Dirichlet boundary conditions, $\alpha =1$, $\beta =0$. For Neumann boundary conditions, $\alpha =0$ and $\beta =1$.

Users need implement the eval() pure virtual function in a derived class. Parameter t is the time variable. Parameter x is the boundary state value either $a$ or $b$. Parameter ip is the index ID of the $i$-th problem when multiple problems are solved together.

Description:

The nested struct PDEParabolic1D::InitialValueFunct is designed to serve as both the initial value function $v(t_0,x)$ in (4) and the update function to support events and resets. We will describe handling events and resets in a separate section.

Users need implement the eval() pure virtual function in a derived class, where

• •

  Parameter x is a $n$-element buffer of discretized spatial states that are ordered from lower boundary $a$ to high boundary $b$.

• •

  Parameter t is the current time $t$. If $t$ is greater than the PDE start time $t_0$, parameter v_x holds the PDE solution $v(t,x)$ computed by a solver up to this current time. If the actual problem requires $v(t,x)$ to be updated at events or resets, v_x should be updated on the return of this function call. The updated $v(t,x)$ then serves as the new initial values for on-going process.

• •

  Parameter event_id indicates if the current time $t$ is on one of the events specified by the user. If so, $t==\text{event\_time[event\_id]}$ in a user supplied event time array event_time. Otherwise, event_id is passed as a negative value.

• •

  Parameter v_x is a $m\times n$ matrix in-out parameter, where $m$ is the number of problems to be solved simultaneously subject to the same PDE, $n$ is the number of spatial discretization states. On call, it holds the PDE solution $v(t,x)$ computed by a solver up to time $t$. On return, it receives updated initial values.

Description:

The nested struct PDEParabolic1D::SourceFunct represents the external source function $s(t,x)$ in (3). Users need implement the add() pure virtual function in a derived class with the following parameters.

• •

  Parameter x is a pointer to the $n$-element buffer of discretized spatial state variable $x$.

• •

  Parameter n specifies the number of spatial states.

• •

  Parameter t is the current time $t$ to evaluate source function $s(t,x)$.

• •

  Parameter ip indicates the source function is to be evaluated for the $i$-th problem when multiple problems governed by the same PDE are solved together.

• •

  Parameter c is a scalar quantity to be multiplied to $s(t,x)$.

• •

  Parameter r is a pointer to a $n$-element vector $𝐫$ where the source function vector $𝐬=s(t,x)$ are scaled by $c$ and then added to $𝐫$, NOT assigned. In other words, vector $𝐫=𝐫+c\times 𝐬$ on function return.

Description:

The nested struct PDEParabolic::Specification groups a list of control setting to specify how the PDE should be solved.

Data member event_t holds important event time points $t_i\in (t_0,T)$ that we want a PDE solver to make a step at $t_i$.

Data member n_problems is the number of problems to be solved together. These problems all follow the same PDE and differ from each other on initial values, boundary conditions and source functions.

Data member t_smooth_initial_values is a time $s$, , where $t_0$ and $T$ are the PDE start and end time respectively. If such a time $s$ is provided, we will solve the problem in two stages. In an initial smoothing stage $t\in (t_0,s)$, a 4-th order accurate finite difference method is used to smooth out initial values from $t_0$ to $s$. A PDE solver then takes $v(s,x)$ as the initial value for the subsequent stage $t\in (s,T)$.

Boolean data member reset_initial_value_at_events is true if we want the solver to reset initial values at those event times. At an event time $t$, the initial value function supplied by a user is called by the PDE solver. Inside the initial value function, users may update the solution $v(t,x)$, which becomes the new initial value to continue the process.

Boolean data member reset_initial_value_at_steps is true if we want the solver to reset initial values at each discretization time step, which effectively treats each step as a new initial value problem. At each time discretization step $t$, user supplied initial value function is called to update $v(t,x)$.

Boolean data member stop_at_last_event instructs a PDE solver to stop at the last event time. Otherwise the solver will run all the way to PDE end time $T$.

Description:

The nested struct PDEParabolic1D::Observation is a helper struct for users and NPL internally to mark observation points on PDE mesh grids, i.e. used in model calibration. Data members t and x are the observation time value $t$ and state variable value $x$. Data member range_v allows users to specify the observed value range.

Description:

The nest struct PDEParabolic1D::Control is an interface for advanced use of PDE solvers. Function $v(t,x)$ could be based on a mathematical model with model parameter $\theta$. Given some observation data $\overline{v}(t_i,x_j)$, we could calibrate model parameter $\theta$ by formulating an optimal control problem, for example,

, where $𝒜(v)$ is the partial derivative operator at the right hand side of (3). In this formulation, PDE constraint (8) translates to a system of nonlinear equality constraints, i.e. a system of difference equations in finite difference methods.

Description:

The rectangular domain $[t_0,T]\times [a,b]$ of function $v(t,x)$ can be queried by the following PDEParabolic1D class member functions:

1. 1.

   Returns the lower state boundary $a$.

2. 2.

   Returns the upper state boundary $a$.

3. 3.

   Returns the PDE start time $t_0$.

4. 4.

   Returns the PDE end time $T$.

Description:

Users need derive from PDEParabolic1D interface class and implement this coefficient() pure virtual member function to evaluate coefficient functions $f(t,x)$, $g(t,x)$ and $h(t,x)$ that appears in (3).

Parameter t is the value of time variable $t$. Parameter x is a pointer to $n$-element buffer of discretized states. Parameter n_x specifies the number $n$ of states in the discretization buffer x.

Parameter _ftx, _gtx and _htx each is a pointer to $n$-element buffer to store $f(t,x)$, $g(t,x)$ and $h(t,x)$ on function return.

For efficiency, an concrete derived class doesn’t have to fill the content of _ftx, _gtx and _htx by zero if the corresponding differential term doesn’t exist in (3). Instead, it could leave the content not touched.

Description:

Create a solver object with given discretization settings. Parameter pde is an object of user derived PDE class under PDEParabolic1D base class. Parameter dx suggests the spatial discretization step size. Parameter dt suggests the time dimension discretization step size.

Description:

1. 1.

   Set the implicitness parameter $0\le \varphi \le 1$ in the general different equation form that is documented in Mathwrist’s white paper. Crank-Nicolson method should set $\varphi =0.5$. $\varphi =0$ or 1 is effectively a fully explicit or a fully implicit method respectively.

2. 2.

   Return the implicitness parameter $\varphi$ used by the solver.

3. 3.

   Return an object reference to the PDE problem to be solved.

4. 4.

   Return a buffer of discretized spatial states that is used by the solver.

5. 5.

   Return the time step size.

Description:

Solve $m$ number of initial value problems simultaneously using Crank Nicolson method.

Parameter f_t is the initial value funtion. It also serves as the callback function to update $v(t,x)$ at events or resets.

Parameter bc_lo and bc_hi specify the lower and upper boundary conditions respectively.

Parameter s_tx is a pointer to the external source function. Pass nullptr if the PDE is homogeneous.

Parameter spec is a problem specification.

Parameter V a $m\times n$ matrix to hold the numerical solutions, where $m$ is the number of problems in specification and $n$ is the number of discretized states.

Description:

Parameter pde is an object of user derived PDE class under PDEParabolic1D base class. Parameter n is the degree of Chebyshev approximation polynomial. Parameter q is the number of error correction iterations used in implicit IDeC, REQUIRED: . Parameter dt suggests the time dimension step size used in implicit IDeC.

Description:

1. 1.

   Return an object reference to the PDE problem to be solved.

2. 2.

   Return the collocation points in physical states. The number of collocation points is equal to the degree of Chebyshev approximation polynomial + 1.

3. 3.

   Return the time step size.

4. 4.

   Return the number of error correction iterations.

Description:

Solve $m$ number of initial value problems simultaneously using spectral collocation method.

Parameter f_t is the initial value funtion. It also serves as the callback function to update $v(t,x)$ at events or resets.

Parameter bc_lo and bc_hi specify the lower and upper boundary conditions respectively.

Parameter s_tx is a pointer to the external source function. Pass nullptr if the PDE is homogeneous.

Parameter spec is a problem specification.

Parameter C is a $m\times (n+1)$ coefficient matrix on return of the function. It holds Chebyshev expansion coefficients to approximate the solutions, where $m$ is the number of problems in specification and $n$ is the degree of Chebyshev approximation polynomial. The numerical solution of each problem can be recovered by creating a ChebyshevCurve object with a row of coefficients of C.

Description:

Enumeration item PDEParabolic_MoL::O2 and PDEParabolic_MoL::O4 instruct a solver to approximate the partial derivative operator $𝒜(v)$ by a 2nd or 4th order accurate finite difference scheme respectively.

Description:

Parameter pde is an object of user derived PDE class under PDEParabolic1D base class. Parameter o is the order of finite difference scheme to be used. Parameter q is the number of error correction iterations used in implicit IDeC, REQUIRED: . Parameter dx suggests the finite difference discretization step size in states. Parameter dt suggests the time dimension step size used in implicit IDeC.

Description:

1. 1.

   Return an object reference to the PDE problem to be solved.

2. 2.

   Return the order of finite difference scheme.

3. 3.

   Return a buffer of discretized spatial states that is used by the solver.

4. 4.

   Return the time step size.

5. 5.

   Return the number of error correction iterations.

Description:

Solve $m$ number of initial value problems simultaneously using method of lines.

Parameter f_t is the initial value funtion. It also serves as the callback function to update $v(t,x)$ at events or resets.

Parameter bc_lo and bc_hi specify the lower and upper boundary conditions respectively.

Parameter s_tx is a pointer to the external source function. Pass nullptr if the PDE is homogeneous.

Parameter spec is a problem specification.

Parameter V a $m\times n$ matrix to hold the numerical solutions, where $m$ is the number of problems in specification and $n$ is the number of discretized states.