Dynamic Programming Method for Best Piecewise Linear Approximation for Vector Field of Nonlinear Boundary Value Problems on the Interval [0, 1]

TL;DR

This paper introduces a dynamic programming approach for optimally approximating vector fields in nonlinear boundary value problems, ensuring boundary condition adherence and numerical stability, with convergence guarantees under finer discretization.

Contribution

It presents a novel dynamic programming method for piecewise linear approximation of vector fields in boundary value problems, emphasizing stability and convergence.

Findings

Method guarantees convergence to true solution with finer discretization

Ensures boundary conditions are maintained during approximation

Provides a stable computational framework for nonlinear BVPs

Abstract

An important problem that arises in many engineering applications is the boundary value problem for ordinary differential equations. There have been many computational methods proposed for dealing with this problem. The convergence of the iterative schemes to a true solution, when one such exists, and their numerical stability are the central issues discussed in the literature. In this paper, we discuss a method for approximating the vector field, maintaining the boundary conditions and numerical stability. If a true solution exists, a subsequence of solutions convergent to one such can be produced, by finer discretization of the solution space.

Equations

$$\left.\begin{array}[]{cl}{\mathbf{x}}^{\prime}(t)={\mathbf{f}}({\mathbf{x}}(t),\,t),&\textrm{ for }0<t<1,\\ \beta({\mathbf{x}}(0),\,{\mathbf{x}}(1))=0,&\end{array}\right\}$$

$$\left.\begin{array}[]{l}{\mathbf{y}}^{\prime}(t)=A_{i}(t){\mathbf{y}}(t)+{\mathbf{b}}_{i}(t),~{}\textrm{ for }t_{i}<t<t_{i+1}\\ \textrm{ with the initial condition }{\mathbf{y}}(t_{i})=\theta_{i}\end{array}\right\}\textrm{ for }0\leq i\leq N-1$$

$$\left.\begin{array}[]{cl}\beta(\theta_{0},\,\theta_{N})~{}=~{}0,&\textrm{ where }\theta_{N}={\mathbf{y}}(1^{-})\\ {\mathbf{y}}(t_{i+1}^{-})~{}=~{}\theta_{i+1},&\textrm{ for }0\leq i\leq N-1,\textrm{ and }\\ A_{i}(t_{i+1})\theta_{i+1}+{\mathbf{b}}_{i}(t_{i+1})&~{}=A_{i+1}(t_{i+1})\theta_{i+1}+{\mathbf{b}}_{i+1}(t_{i+1}),\\ &\textrm{ for }0\leq i\leq N-2\end{array}\right\}$$

$${\mathscr{\LARGE E}}=\sum_{i=0}^{N-1}\int_{t_{i}}^{t_{i+1}}\mathcal{G}{\bigg{(}}A_{i}(t){\mathbf{y}}(t)+{\mathbf{b}}_{i}(t)-{\mathbf{f}}({\mathbf{y}}(t),\,t){\bigg{)}}\,\mathrm{d}t$$

$${\mathbf{x}}^{\prime }(t)=\mathbf{f}(\mathbf{x}(t),t),\text{ for }0<t<1,$$

$${\beta }_0(\mathbf{x}(0))={\beta }_1(\mathbf{x}(1))=0,$$

$$\left.\begin{array}[]{cl}\beta_{0}(\theta_{0})\,=\,0,~{}\beta_{1}(\theta_{N})\,=\,0,&\textrm{ where }\theta_{N}={\mathbf{y}}(1)\\ {\mathbf{y}}(t_{i+1}^{-})~{}=~{}\theta_{i+1},&\textrm{ for }0\leq i\leq N-1,\textrm{ and }\\ A_{i}(t_{i+1})\theta_{i+1}+{\mathbf{b}}_{i}(t_{i+1})&~{}=~{}A_{i+1}(t_{i+1})\theta_{i+1}+{\mathbf{b}}_{i+1}(t_{i+1}),\\ &\textrm{ for }0\leq i\leq N-2\end{array}\right\}$$

$${\mathbf{y}}(t)=\textit{\Large e}^{{\bigg{[}}\int_{t_{i}}^{t}A_{i}(s)\,\,\mathrm{d}s{\bigg{]}}}\,\theta_{i}+\int_{t_{i}}^{t}\textit{\Large e}^{{\bigg{[}}\int_{s}^{t}A_{i}(v)\,\,\mathrm{d}v{\bigg{]}}}\,{\mathbf{b}}_{i}(s)\,\,\mathrm{d}s,$$

$$\displaystyle{\Large{\varepsilon}}_{k}~{}=~{}\int_{t_{k-1}}^{t_{k}}\mathcal{G}{\bigg{(}}A_{k-1}(t){\mathbf{y}}(t)+{\mathbf{b}}_{k-1}(t)-{\mathbf{f}}({\mathbf{y}}(t),\,t){\bigg{)}}\,\mathrm{d}t$$

$${\mathcal{E}}_k={\mathcal{E}}_{k-1}+{\epsilon }_k$$

$$\left.\begin{array}[]{c}\beta_{0}(\theta_{0})\,=\,0,~{}\textrm{ if }k=1\\ \beta_{N}(\theta_{N})\,=\,0,~{}\textrm{ if }k=N\\ {\mathbf{y}}(t_{k}^{-})~{}=~{}\theta_{k},~{}\textrm{ and }\\ A_{k-1}(t_{k})\theta_{k}+{\mathbf{b}}_{k-1}(t_{k})\,=\,A_{k}(t_{k})\theta_{k}+{\mathbf{b}}_{k}(t_{k}),~{}\textrm{ if }k<N\\ \end{array}\right\}$$

$${\Gamma }_k=\{({\pi }_0,{\pi }_1,\dots ,{\pi }_k)|({\pi }_{i-1},{\pi }_i)\in {\Omega }_i,\text{ for }1\le i\le k\}$$

$$S_{k}(\pi_{k})\hskip 14.22636pt=\min_{\scriptsize{\begin{array}[]{c}(\pi_{0},\ldots,\pi_{k-1})\\ ~{}\textrm{\footnotesize such that}\\ (\pi_{0},\ldots,\pi_{k-1},\pi_{k})\,\in\,\Gamma_{k}\end{array}}}{\mathscr{\LARGE E}}_{k}(\pi_{0},\ldots,\pi_{k-1})$$

$$\underset{({\pi }_0,\dots ,{\pi }_{k-1},{\pi }_k)\in {\Gamma }_k}{\min }{\mathcal{E}}_k({\pi }_0,\dots ,{\pi }_{k-1})=\underset{{}_{\text{and }({\pi }_{k-1},{\pi }_k)\in {\Omega }_k}^{({\pi }_0,\dots ,{\pi }_{k-1})\in {\Gamma }_{k-1}}}{\min }\{{\mathcal{E}}_{k-1}({\pi }_0,\dots ,{\pi }_{k-2})+{\epsilon }_k({\pi }_{k-1})\}$$

$$S_k({\pi }_k)=\underset{{}_{({\pi }_{k-1},{\pi }_k)\in {\Omega }_k}^{{\pi }_{k-1}\text{such that}}}{\min }\{S_{k-1}({\pi }_{k-1})+{\epsilon }_k({\pi }_{k-1})\}$$

$${\pi }_0^{\ast }({\pi }_1)=\{{\pi }_0|({\pi }_0,{\pi }_1)\in {\Omega }_1\text{ and }\epsilon ({\pi }_0)=S_1({\pi }_1)\}$$

$${\pi }_{k-1}^{\ast }({\pi }_k)=\{{\pi }_{k-1}|({\pi }_{k-1},{\pi }_k)\in {\Omega }_k\text{ and }{S}_k({\pi }_k)=S_{k-1}({\pi }_{k-1})+{\epsilon }_k({\pi }_{k-1})\}$$

$$\left.\begin{array}[]{c}\exists{\xi_{N}}:\,\beta(\theta_{0},\xi_{N})\,=\,0,~{}\textrm{ if }k=1\\ \exists{\xi_{0}}:\,\beta(\xi_{0},\theta_{N})\,=\,0,~{}\textrm{ if }k=N\\ {\mathbf{y}}(t_{k}^{-})~{}=~{}\theta_{k},~{}\textrm{ and }\\ A_{k-1}(t_{k})\theta_{k}+{\mathbf{b}}_{k-1}(t_{k})\,=\,A_{k}(t_{k})\theta_{k}+{\mathbf{b}}_{k}(t_{k}),~{}\textrm{ if }k<N\\ \end{array}\right\}$$

$$S_k({\pi }_0,{\pi }_k)=\underset{{}_{({\pi }_0,{\pi }_1,\dots ,{\pi }_{k-1},{\pi }_k)\in {\Gamma }_k}^{({\pi }_1,\dots ,{\pi }_{k-1})\text{such that }}}{\min }{\mathcal{E}}_k({\pi }_0,{\pi }_1,\dots ,{\pi }_{k-1})$$

$$S_k({\pi }_0,{\pi }_k)=\underset{{}_{({\pi }_{k-1},{\pi }_k)\in {\Omega }_k}^{{\pi }_{k-1}\text{such that}}}{\min }\{S_{k-1}({\pi }_0,{\pi }_{k-1})+{\epsilon }_k({\pi }_{k-1})\}$$

$$\min_{\scriptsize{\begin{array}[]{c}(\pi_{0},\,\pi_{N})\\ \textrm{ such that }\\ \beta(\theta_{0},\,\theta_{N})=0\\ \end{array}}}S_{N}(\pi_{0},\pi_{N})~{}~{}=\min_{\scriptsize{\begin{array}[]{c}(\pi_{1},\,\ldots,,\,\pi_{N-1})\\ \textrm{ such that }\\ (\pi_{0},\,\pi_{1},\,\ldots,\,\pi_{N})\,\in\,{\mathcal{V}}_{N}\bigcap\Gamma_{N}\\ \end{array}}}{\mathscr{\LARGE E}}_{N}(\pi_{0},\pi_{1},\ldots,\pi_{N-1})$$

$${\pi }_{k-1}^{\ast }({\pi }_0,{\pi }_k)$$

$$\begin{array}[]{rccl}(\pi_{0}^{\ast},\,\pi_{N}^{\ast})&=&\textrm{argmin}&S_{N}(\pi_{0},\,\pi_{N})\\ &&\scriptstyle{(\pi_{0},\,\pi_{N})}~{}\textrm{\footnotesize such that}&\\ &&\scriptstyle{B(\theta_{0},\,\theta_{N})\,=\,0}&\\ \end{array}$$

$${\mathbf{x}}_0^{\prime }(t)\approx \mathbf{f}({\mathbf{x}}_0(t),t),\text{ for }0<t<1$$

$$\left.\begin{array}[]{c}{\mathbf{x}}_{k}^{\prime}(t)={\mathbf{f}}_{k}(t),~{}\textrm{ for }0<t<1\\ \textrm{with the initial condition }{\mathbf{x}}_{k}(0)=\theta_{0}\end{array}\right\}$$

$$\left.\begin{array}[]{c}{\mathbf{x}}_{k+1}^{\prime}(t)={\mathbf{f}}({\mathbf{x}}_{k+1}(t),\,t)~{}\textrm{ for }0<t<1\\ \textrm{with the initial condition }{\mathbf{x}}_{k+1}(0)=\theta_{0}\end{array}\right\}$$

$$\left.\begin{array}[]{c}{\mathbf{y}}_{k}^{\prime}(t)={\mathbf{f}}({\mathbf{x}}_{k+1}(t),\,t)-{\mathbf{f}}_{k}(t),~{}\textrm{ for }0<t<1\\ \textrm{with the initial condition }{\mathbf{y}}_{k}(0)=0\end{array}\right\}$$

$$\mathbf{f}({\mathbf{x}}_{k+1}(t),t)$$

$$\left.\begin{array}[]{c}{\mathbf{y}}_{k}^{\prime}(t)=A_{k}(t){\mathbf{y}}_{k}(t)+{\mathbf{b}}_{k}(t),~{}\textrm{ for }0<t<1\\ \textrm{with the initial condition }{\mathbf{y}}_{k}(0)=0\end{array}\right\}$$

$${\mathbf{y}}_{k}(t)=\int_{0}^{t}\textit{\Large e}^{{\bigg{[}}\int_{s}^{t}A_{k}(v)\,\,\mathrm{d}v{\bigg{]}}}\,{\mathbf{b}}_{k}(s)\,\,\mathrm{d}s$$

$${\mathbf{y}}_k(t)=-A_k^{-1}(t){\mathbf{b}}_k(t),0\le t\le 1$$

Full text

1 Introduction

During the last few decades there has been a remarkable growth of interest in problems associated with solving linear and nonlinear ordinary differential equations satisfying boundary conditions. For many of the nonlinear boundary value problems that occur in engineering and applied sciences, it is difficult to obtain a solution analytically. For a nonlinear boundary value problem, the difficulty lies in establishing the existence of a solution mathematically, though in some cases multiple solutions exist. Approximation of the solution space of a given differential equation has gained importance as it speeds up or helps in solving the problem efficiently. These approximation methods can be put into two classes: (i) those in which a solution is approximated numerically at a number of discrete points of the domain, and (ii) those in which a solution is approximated by a finite number of terms of a sequence of functions. The approach in (ii) is called a weighted residual method. In most of the numerical methods described in the literature, we may have to add a sufficient number of some undetermined variables with implicitly assumed conditions at one end of the domain, and adjust the additional variables until the required conditions are satisfied at the other end to obtain the solution of the boundary value problem [[2] – [4], [7], [10] – [13]], and further approximate the derivatives of the dependent variables with forward, backward or central difference operators defined at the grid points [9]. In the first part of a numerical scheme as just mentioned, the convergence may be very slow, and in the second part, convergence and stability of the particular difference scheme may depend on the selection of the approximations used for the derivatives involved in the differential equation and the boundary conditions. It may also be the case that the chosen difference method is not numerically stable, resulting in chaos phenomenon creeping into the iterative schemes, at places where matrix inversions are utilized, without the solver being explicitly aware of its entry. In some cases, suitable regularization and relaxation conditions, involving more variables, may have to be added to the constraints formulated in the previous steps. In the weighted residual methods, difference equations are generated using approximation methods with piecewise polynomial solutions [9].

Among the most popular and successful techniques for solving boundary value problems with nonlinearities is Galerkin procedure. In this approach, the solution of the ordinary differential equation is expressed as a linear combination of certain basis functions, and the coefficients of the basis functions are determined by requiring that the residual be orthogonal to each of the basis functions. The difficulty lies in the selection of basis functions to obtain the desired solution, that can take care of the boundary conditions simultaneously. In recent times, the concept of piecewise linear approximation of the differential equation gained momentum [[5], [7], [9] – [10]] The two point boundary value problems are approximated by piecewise linear ones which have analytical solutions and reduced to finding the slope of the solution at the left boundary, so that the boundary conditions at the right end of the interval are satisfied. This approach results in a complex system of non-linear algebraic equations. Some more recent and highly efficient algorithms [6], for solving these complex systems of differential and algebraic equations (DAE) can be used for computational purpose.

The motivation for the present work is the necessity of a new method that is applicable for most or all general continuous vector fields and general boundary conditions. The objective of our study is to find efficiently a solution, if exists, by an algorithm, such that the algorithm is able to detect and report when there is no solution, if an error term does not fall below a threshold, despite using various approximation schemes with several basis functions. We propose a new method based on dynamic programming for solving boundary value problems in one variable. The dynamic programming based formulation is adapted for obtaining an optimal approximation for the vector field for genera of two-point boundary value problem, which is usually formulated as an optimal control problem in the literature [1]. For improving an initial approximation, repeated application of the dynamic programming algorithm with refined discretization of the parameter space can be used. A modified Newton-Raphson method for improving an approximate solution along with updating the initial value is also discussed. These aspects are newly introduced in this work. Possible extensions to formulation of solution methods for optimal control problem and for boundary value problems of partial differential equations defined on compact simply connected convex domains are briefly described.

2 Best Piecewise Linear Approximations for Vector Fields of Two-point Boundary Value Problems

In this section, we consider the boundary value problem

[TABLE]

where ${\mathbf{f}}\in{\mathscr{\LARGE C}}({\mathbb{R}}^{n}\times[0,\,1];\,{\mathbb{R}}^{n})$, $\beta\in{\mathscr{\LARGE C}}({\mathbb{R}}^{n}\times{\mathbb{R}}^{n};\,{\mathbb{R}})$, and ${\mathbf{x}}^{\prime}(t)=\frac{\,\mathrm{d}{\mathbf{x}}}{\,\mathrm{d}t}$. The vectors are written as columns of appropriate dimensions. The objective of obtaining a piecewise linear approximation is stated as follows:

Let, for some fixed integer $N\geq 2$, $0=t_{0}<t_{1}<\ldots<t_{N}=1$ be fixed node points of the time interval $[0,\,1]$. 2. 2.

Then, we obtain optimal values for the parameters of form $\pi_{i}=(A_{i}(t),\,{\mathbf{b}}_{i}(t),\,\theta_{i})$, $0\leq i\leq N-1$, where $A_{i}$ and ${\mathbf{b}}_{i}$ are $n\times n$ and $n\times 1$ matrices consisting of undetermined linear combination of some fixed basis functions defined on $[t_{i},\,t_{i+1}]$, and $\theta_{i}$ is $n\times 1$ column vector of undetermined constants, such that the curve ${\mathbf{y}}(t)$ satisfying

[TABLE]

subject to the boundary and continuity conditions

[TABLE]

minimizes the error defined by

[TABLE]

where the integrand $\mathcal{G}$ can be chosen to be any nonnegative continuous function defined on ${\mathbb{R}}^{n}$ and equal to 0 at the origin.

We first consider a special case of (2.1), in which the boundary conditions are variables seperable (i.e. of the form $\beta_{0}({\mathbf{x}}(0))=0$ and $\beta_{1}({\mathbf{x}}(1))=0$, where $\beta_{0},\,\beta_{1}\in{\mathscr{\LARGE C}}({\mathbb{R}}^{n};\,{\mathbb{R}})$). Later, we shall indicate the modifications needed to extend the method to the general case of boundary conditions described in (2.1).

2.1 Piecewise Linear Approximation for the Special Case

In this section, we consider the boundary value problem

[TABLE]

where ${\mathbf{f}}\in{\mathscr{\LARGE C}}({\mathbb{R}}^{n}\times[0,\,1];\,{\mathbb{R}}^{n})$, and $\beta_{0},\,\beta_{1}\in{\mathscr{\LARGE C}}({\mathbb{R}}^{n}\times{\mathbb{R}}^{n};\,{\mathbb{R}})$, where ${\mathscr{\LARGE C}}(X;\,Y)$ is the space of continuous functions from $X$ into $Y$ and ${\mathbb{R}}$ is the set of real numbers. The objective being described in the sequel is to find optimal values for the parameters $\pi_{i}=(A_{i},\,B_{i},\,\theta_{i})$, $0\leq i\leq N-1$, for which the error functional in (2.3) is minimized by the solution satisfying (2) subject to the boundary and continuity conditions

[TABLE]

The solution for (2) in the interval $[t_{i},\,t_{i+1})$, $0\leq i\leq N-1$, is

[TABLE]

which is expressible explicitly in terms of the parameter values. Substituting the expression for ${\mathbf{y}}(t)$ in (2.3), the value of ${\mathscr{\LARGE E}}$ can be found for any prescribed values of the parameters. If the parameter space is discretized, the values of ${\mathscr{\LARGE E}}$ can be tabulated for various values of the parameters. However, using the additive property of the error functional ${\mathscr{\LARGE E}}$, it is possible to formulate a dynamic programming method for obtaining the tables quickly and efficiently. For this purpose, we define a $k$-step error functional ${\mathscr{\LARGE E}}_{k}$, for $0\leq k\leq N$, by the relation ${\mathscr{\LARGE E}}_{0}=0$ and for $k=1,2,\ldots,N$,

[TABLE]

Thus ${\mathscr{\LARGE E}}_{N}$ in (2.8) is actually ${\mathscr{\LARGE E}}$ in (2.3). The functional ${\mathscr{\LARGE E}}_{k}$ depends only on the parameters $(\pi_{0},\pi_{1},\ldots,\pi_{k-1})$. To make this dependency explicit, we sometimes write ${\mathscr{\LARGE E}}_{k}(\pi_{0},\pi_{1},\ldots,\pi_{k-1})$ for ${\mathscr{\LARGE E}}_{k}$; similarly, ${\Large{\varepsilon}}_{k}$ is sometimes written as ${\Large{\varepsilon}}_{k}(\pi_{k-1})$. Now let for $1\leq k\leq N$, $\Omega_{k}$ be the set of parameters $(\pi_{k-1},\,\pi_{k})$ satisfying

[TABLE]

where if $k=N$, the parameter $\pi_{N}$ is interpreted as $\theta_{N}$, and ${\mathbf{y}}(t)$ in the interval $[t_{k-1},\,t_{k})$ is given by (2.6) with $i=k-1$. Let

[TABLE]

Now we define a finite sequence of functions on the parameter space as follows: let $S_{0}(\pi_{0})=0$, and for $1\leq k\leq N$,

[TABLE]

Then, observing that for each fixed $\pi_{k}$,

[TABLE]

we can conclude that

[TABLE]

The value of $\min_{\pi_{N}}S_{N}(\pi_{N})=\min_{\theta_{N}}S_{N}(\theta_{N})$ gives the minimum value of ${\mathscr{\LARGE E}}_{N}$ over the parameter space constrained by the boundary and continuity conditions (2.5).

Now we describe a tabulation procedure for computation of optimal parameters as follows:

Initially compute $S_{1}(\pi_{1})=\min_{(\pi_{0},\,\pi_{1})\in\Omega_{1}}{\Large{\varepsilon}}_{1}(\pi_{0})$, and set

[TABLE] 2. 2.

Now for $2\leq k\leq N$, compute $S_{k}(\pi_{k})=\min_{(\pi_{k-1},\,\pi_{k})\in\Omega_{k}}S_{k-1}(\pi_{k-1})+{\Large{\varepsilon}}_{k}(\pi_{k-1})$, and set

[TABLE]

At the end of the procedure, we obtain a table for $S_{N}(\pi_{N})$ for various values of $\pi_{N}$ (which is $\theta_{N}$). The optimal parameters are found by backtracking as follows:

First find $\pi_{N}^{\ast}=\textrm{argmin}_{\pi_{N}}S_{N}(\pi_{N})$. 2. 2.

Then for $k=N-1,\,\ldots,\,1,\,0$, choose $\pi_{k}^{\ast}$ as an element of $\pi_{k}^{\ast}(\pi_{k+1}^{\ast})$.

The tuple $(\pi_{0}^{\ast},\,\pi_{1}^{\ast},\,\ldots,\,\pi_{N}^{\ast})$ thus obtained are optimal parameters with respect to the error functional (2.3) subject to the boundary and continuity conditions in (2.5).

2.2 Piecewise Linear Approximation for the General Case

Now we describe a generalization of the method of the previous section to the general case of boundary conditions in (2.2). Let for $1\leq k\leq N$, $\Omega_{k}$ be the set of parameters $(\pi_{k-1},\,\pi_{k})$ satisfying

[TABLE]

Let $\Gamma_{k}$ be as in (2.10), but with $\Omega_{k}$ as in (2.13). Now the parameter space constrained by (2.2) is given by $\Gamma_{N}\bigcap{\mathcal{V}}_{N}$, where ${\mathcal{V}}_{N}=\{(\pi_{0},\,\pi_{1},\,\ldots,\,\pi_{N})\,|\;\beta(\theta_{0},\,\theta_{N})=0\}$. However, $S_{k}$ is now defined on $(\pi_{0},\,\pi_{k})$:

[TABLE]

so that

[TABLE]

for $2\leq k\leq N$. Now

[TABLE]

The computation of the forward tables is as follows:

Initially compute $S_{1}(\pi_{0},\,\pi_{1})$ for various values of $(\pi_{0},\,\pi_{1})\in\Omega_{1}$. 2. 2.

Now for $2\leq k\leq N$, compute $S_{k}(\pi_{0},\,\pi_{k})$ from (2.15), and set

[TABLE]

From the table thus computed, the optimal values of the parameters are extracted by backtracking as follows:

First find $(\pi_{0}^{\ast},\,\pi_{N}^{\ast})$ from

[TABLE] 2. 2.

Then for $k=N-1,\,\ldots,\,2,\,1$, choose $\pi_{k}^{\ast}$ as an element of $\pi_{k}^{*}(\pi_{0}^{\ast},\,\pi_{k+1}^{\ast})$.

The tuple $(\pi_{0}^{\ast},\,\pi_{1}^{\ast},\,\ldots,\,\pi_{N}^{\ast})$ thus obtained minimizes the error functional (2.3) subject to the boundary and continuity conditions in (2.2).

3 Improving the Initial Approximation

In this section, we describe a method similar to Newton-Raphson method for improving an initial approximation for a BVP or in particular, for an initial value problem. The correction of the approximation consists of two parts: in the first part, a correction term with 0 initial value for improving the solution in the interior is obtained, and the second part finds optimal correction term in the initial value so that the updated solution satisfies the boundary conditions more accurately. Alternatively, it is also possible to achieve the same by either Picard’s successive approximation (with suitable correction term in the initial value) or repeated application of the dynamic programming method described in the last section with finer discretization of the parameter space restricted to a tube-like set around the initial approximate solution.

3.1 Improving an Approximate Solution with Intial Value Fixed

We assume that the vector field ${\mathbf{f}}$ is differentiable with continuous derivatives of upto as high an order (upto two) as necessary. Let ${\mathbf{x}}_{0}(t)$ be an initial approximate solution satisfying

[TABLE]

The objective is to formulate an efficient method for improving the approximation in (3.1) by successive iterations. Let for $k=0,\,1,\,2,\,\ldots$,

[TABLE]

At step $k+1$, the correction term ${\mathbf{y}}_{k}={\mathbf{x}}_{k+1}-{\mathbf{x}}_{k}$ must be found such that

[TABLE]

Subtracting (3.2) from (3.3), we find that ${\mathbf{y}}_{k}$ must satisfy

[TABLE]

Now using the Taylor series approximation with respect to the first variable (i.e. ${\mathbf{x}}$) upto the first order for ${\mathbf{f}}({\mathbf{x}}_{k+1}(t),\,t)$ we find

[TABLE]

Letting ${\mathbf{b}}_{k}(t)={\mathbf{f}}({\mathbf{x}}_{k}(t),\,t)-{\mathbf{f}}_{k}(t)$ and using (3.1) in (3.4), an approximate correction term is found by solving

[TABLE]

The solution of (3.6) is given by

[TABLE]

The iteration converges fast (at almost quadratic rate) provided the initial approximate solution (3.1) is sufficiently close to the exact solution. If the matrix $A_{k}(t)$ is uniformly boundedly invertible in the interval $[0,\,1]$, then ${\mathbf{y}}_{k}(t)$ can also be taken as

[TABLE]

which is the well-known method for solving for ${\mathbf{x}}$ from ${\mathbf{f}}({\mathbf{x}},t)=0$ for each $0\leq t\leq 1$. This method also modifies the initial condition. If the initial value is required to be updated independently, then we have to find ${\mathbf{y}}_{k}(t)$ from (3.6), with solution given by (3.7). It may be observed when either of (3.7) and (3.8) converges, $\lim_{k\rightarrow\infty}{\mathbf{b}}_{k}(t)=0$, almost everywhere (a.e.), for $0\leq t\leq 1$, which implies, if convergent, either iteration leads to the final solution satisfying ${\mathbf{x}}^{\prime}(t)={\mathbf{f}}({\mathbf{x}}(t),\,t)$ a.e., for $0\leq t\leq 1$.

3.2 Improving the Initial Value

Suppose that the initial value of approximate correction term ${\mathbf{y}}_{k}$ is $\eta_{k}$ (to be determined). Then the correction term is given by

[TABLE]

where $A_{k}(t)$ is as in (3.1) and ${\mathbf{b}}_{k}(t)$ in the following line. Then the parameter $\eta_{k}$ is found based on an optimality criterion. The objective can be formulated as follows:

[TABLE]

If the initial approximation is sufficiently accurate at $t=0$, then $\eta_{k}$ must be small. Thus we can expect that the minimum in (3.10) is attained for $\eta_{k}\approx 0$. The update value can be chosen to be proportional to the gradient of $\mathscr{F}(\eta_{k})$ at 0. Specifically, we can choose $\eta_{k}^{\ast}=-h_{k}\cdot\nabla\mathscr{F}(0)$, for a small positive number $h_{k}$, resulting in an easy updation of the initial value. The multiplier $h_{k}$ can be found by binary search method over an interval of the form $[0,\,M]$, for some sufficiently large positive constant $M$. Alternately it is also possible to find $\eta_{k}^{\ast}$ by solving $\nabla{\mathscr{F}}(\eta)=0$, which gives the iterative formula $\eta_{0}=0$ and for $k=1,\,2,\,\ldots$,

[TABLE]

where $H(\eta)$ is the Hessian matrix of ${\mathscr{F}}(\eta)$. The method is fast and does not require a seperate search for the multiplier constant as in the case of gradient descent method.

3.3 Combining Both: Improving an Initial Approximation by Successive Iterations

In this section we briefly describe a gradient descent method for imporving an initial solution of the boundary value problem (2.1). For this purpose, we assume that $\beta$ is continuously differentiable having continuous derivatives upto Hessian. The update in approximation ${\mathbf{x}}_{k}(t)$ of the $k$-th iteration is given by (3.9) with initial condition ${\mathbf{y}}_{k}(0)=\eta_{k}$, and ${\mathbf{y}}^{\prime}_{k}(t)$ satisfying (3.4). The initial value $\eta_{k}$ of ${\mathbf{y}}_{k}(t)$ is determined by solving the boundary condition. Specifically let $\theta_{0}^{(k)}={\mathbf{x}}_{k}(0)$ and $\theta_{N}^{(k)}={\mathbf{x}}_{k}(1)$. The initial value $\theta_{0}^{(k+1)}={\mathbf{x}}_{k+1}(0)$ is $\theta_{0}^{(k)}+\eta_{k}$, and the final value $\theta_{N}^{(k+1)}$ is ${\mathbf{x}}_{k}(1)+{\mathbf{y}}_{k}(1)$. The updated values $(\theta_{0}^{(k+1)},\,\theta_{N}^{(k+1)})$ must be found such that

[TABLE]

Now

[TABLE]

Substituting the values of $\theta_{0}^{(k+1)}$ and $\theta_{N}^{(k+1)}$ from (3.13) and (3.14) into (3.12), we find

[TABLE]

As (3.15) is to be solved for the vector $\eta_{k}$ from only one equation, we propose first a gradient descent method for minimization of $\beta^{2}$. Further, if $\eta_{k}$ is small, we can evaluate the gradient of $\beta$ in (3.15) with respect to $\eta_{k}$ for $\eta_{k}=0$. Thus $\eta_{k}$ is chosen such that

[TABLE]

where $\beta(\theta_{0},\,\theta_{N})$ and its partial derivatives $\partial_{\theta_{0}}\beta(\theta_{0},\,\theta_{N})$ and $\partial_{\theta_{N}}\beta(\theta_{0},\,\theta_{N})$ are evaluated for $\theta_{0}=\theta_{0}^{(k)}$ and \theta_{N}=\theta_{N}^{(k)}+\int_{0}^{1}\textit{\Large e}^{{\bigg{[}}\int_{s}^{1}A_{k}(v)\,\,\mathrm{d}v{\bigg{]}}}\cdot\,{\mathbf{b}}_{k}(s)\,\,\mathrm{d}s, and $h_{k}\geq 0$ is a small multiplier that can be found, for example, by binary search method in the interval $[0,\,M]$ for some constant $M>0$. However, it is important to constrain $h_{k}$ to be close to [math], since $\eta_{k}$ must be restricted such that ${\mathbf{y}}_{k}(t)$ never leaves a tube-like set that can be determined for convergence of the Newton-Raphson method. Alternately, it is also possible to find $\eta^{\ast}$ such that $\nabla_{\eta}\beta^{2}=0$, where $\beta$ is as in (3.15). In this case, the update in $\eta_{k}$ is found by the following iteration: $\eta_{0}=0$ and for $k=1,\,2,\,\ldots$,

[TABLE]

where $H(\eta)$ is the Hessian matrix of $\beta^{2}$ with $\beta$ as in (3.15), $\theta_{0}=\theta_{0}^{(k)}$ and $\theta_{N}=\theta_{N}^{(k)}$. The method is fast and converges to the true boundary values provided the initial approximation is sufficiently accurate.

4 Summary of Boundary Value Problem and its Extension to Optimal Control Problem and Multidimensional Cases

In this section, we summarize the dynamic programming method, with an error function using uniform metric. The dynamic programming method can be extended to a vary large class of metrics. We bring out the essential characteristic that is needed for formulation of a dynamic programming based discrete optimization method. A special formulation for the optimal control problem is given.

Towards the end of the section, we briefly indicate how to extend the dynamic programming method for the partial differential equations defined in the interior of a simply connected compact, preferably convex, domain, with a regular boundary together with prescribed conditions on it that a solution must satisfy. It is assumed that the boundary is prescribed by a covering as in an atlas. Applications of this method for solving partial differential equations can be found in remote sensing, spectroscopy and tomography, in which regions of physical matter of different permeating, penetrating, reflexivity or resistivity properties, that can affect a flow field, are estimated by the modeling parameters, using measurements taken at the surface and some interior points, where the measurements at these interior points may be assumed or default values, and comparing the reconstruction with another model, which is assumed to be free from anomalies.

4.1 Summary of Dynamic Programming Method for the Solution to a Boundary Value Problem

The method described in Section 2 can be applied with any nonnegative error functional ${\mathscr{\LARGE E}}$ for which the $k$-step error functional ${\mathscr{\LARGE E}}_{k}$ can be evaluated based on ${\mathscr{\LARGE E}}_{k-1}$ and ${\Large{\varepsilon}}_{k}$. In particular, if for some function $g$ defined on $[0,\,\infty)\times[0,\,\infty)$, the following holds

[TABLE]

and ${\mathscr{\LARGE E}}_{0}$ is specified, then the dynamic programming method can be still applied with the function $g$ in stead of $+$.

One of the important error functionals is with respect to uniform metric, which takes the form

[TABLE]

The $k$-step error funcitonal for (4.2) is given by

[TABLE]

where $1\leq k\leq N$. The error functional (4.2) is best suited for obtaining an initial approximation, followed by the Newton-Raphson method of iterative improvement. A bound on the error in the initial approximation for convergence of the Newton-Raphson mothod for computing the solution of algebraic or analytic equations can be explicitly found. Now as the formulation for the differential equation that the error correction term ${\mathbf{y}}_{k}$ satisfies in Sections 3.1 and 3.2 is obtained by treating the variable $t$ fixed, the same bound works for convergence of ${\mathbf{y}}_{k}$ (or ${\mathbf{y}}_{k}^{\prime}$).

Further, with the error functional given by (4.2) the dynamic programming method can be applied repeatedly until the desired precision is achived using finer and finer discretizations and restricting the search to only the tube-like set around the previously obtained optimal parameters. Decisions concerning which parameters (not necessarily optimal with respect to $S_{k}$) to retain in subsequent iterations can be made based on quantitative measures such as stiffness at the parameter value for identifying the tube set locally. A measure similar to stiffness is the spectrum of the matrix $A_{k}(t)$ defined in (3.1), which indicates as the iteration progresses which components move inward and which outward of the tube. The spectrum of $A_{k}^{\prime}(t)$ indicates the torsion and oscillatory properties of the solution. These measures can be evaluated at the parameter value from the vector field ${\mathbf{f}}$ without requiring complete or even part of the actual solution. Thus besides the tables for $S_{k}$, auxiliary tables containing information regarding stiffness or other measures can be used for choosing the parameter values in the dynamic programming algorithm. The usefulness of the auxiliary tables is especially significant when using larger discretization step-sizes. The auxiliary information can expedite the search method by eliminating unwanted parameter values and retaining only those parameter values that could actually produce the true solution, so that in the subsequent iterations, closer approximations (i.e., with smaller error) are produced. The basic abstract model for the error functional in (4.1) can be recast in such a way that allowances for parameter dependences in the accumulating function $g$ and dependence of the $k$-step cost function ${\mathscr{\LARGE E}}_{k}$ on a future state to reach, which is described by the model parameter $\pi_{k}$, for $1\leq k\leq N$, are explicated, as follows:

[TABLE]

and ${\mathscr{\LARGE E}}_{0}(\pi_{0})$ is specified. The model parameter sequence $(\pi_{0},\,\ldots,\,\pi_{k-1},\,\pi_{k})$ in the arguments of the $g$-function in (4.4) allows the designer to take into consideration costs incurred due to lag or drag involved in the course from an initial or past state up to until the current state, described by the model parameter $\pi_{k-1}$, to reach the next state described by the model parameter $\pi_{k}$, for $1\leq k\leq N$.

4.2 Extension to the Solution of Optimal Control Problem

In this section, we describe a method to solve the optimal control problem

[TABLE]

where ${\mathbf{f}}\in{\mathscr{\LARGE C}}({\mathbb{R}}^{n}\times[0,\,1];\,{\mathbb{R}}^{n})$, $\beta\in{\mathscr{\LARGE C}}({\mathbb{R}}^{n}\times{\mathbb{R}}^{n};\,{\mathbb{R}})$, and ${\mathcal{E}}({\mathbf{u}})$ is a cost functional that can be written as sum of cost functionals over the intervals $[t_{i},\,t_{i+1}]$, $0\leq i\leq N-1$. Taking ${\mathbf{u}}(t)=A_{i}(t){\mathbf{x}}(t)+b_{i}(t)-{\mathbf{f}}({\mathbf{x}}(t),\,t)$, on $[t_{i},\,t_{i+1}]$, $0\leq i\leq N-1$, the dynamic programming method and improvements of approximations described in the previous section can be utilized to find an approximate optimal solution for the control function ${\mathbf{u}}(t)$. For other applications, for instance, in [14], the optimal control based formulation is used for estimation of sets reachable from an initial state.

4.3 Extension for Boundary Value Problems for Partial Differential Equations Defined over a Simply Connected Compact Domain of the Euclidean Space

In this subsection, we describe a method to solve a boundary value problem of a partial differential equation, defined over a simply connected compact domain, which is a subset of ${\mathbb{R}}^{d}$, for some positive integer $d\geq 2$. For simplicity of description, the domain is further assumed to be convex, although this assumption is not always necessary. The boundary is assumed to be sufficiently regular and specified by smooth surfaces, parameterized as in an atlas, such as in the case of a sphere. The boundary is then propagated inward by, for example, computing the Euclidean distance from a point inside the domain. This approach is called boundary propagation or front propagation [8]. A change of coordinates, that is consistent with the description of the domain as consisting of concentric surfaces diffeomorphic to the initial boundary of the domain, is performed in the given partial differential equation. These surfaces are parametrically described in such a way that some geometric attribute remains constant for each surface, and hence they are called level sets of the propagating boundary. If the initial boundary is convex, with its unit normal pointing outward, then the boundary can be propagated inward by collecting points obtained by subtracting vectors normal to the boundary at the initial points on the boundary with small multipliers as absolute value from those points on the initial boundary, as in the gradient descent method, with varying multipliers. But this propagation may not eventually end in a single point, and, depending on the absolute values chosen as multipliers of the surface normal vectors, various shapes can be realized. In order to construct a set of concentric surfaces ending in a point, the Euclidean distance is computed from a point, which is termed as the center of the domain, to various points on the boundary, whose normal always points outwards, by assumption. A criterion for an interior point to be the center is stated as follows: let the distance of an interior point to the boundary be defined to be the maximum distance from it to any point on the boundary, and let the center be chosen to be the interior point that results in a minimum distance to the boundary among all interior points. This method of choosing the center is appropriate with compact convex regular boundaries with simply connected interior domains. Then, from the center, the distance to a point on the boundary is multiplied by a scale, as in the case of projective coordinate system, but the scale parameter is chosen to be a real number between 0 and 1. The level sets are the surfaces corresponding to the same scale. Now, appropriate dynamic programming tables are constructed that approximate a solution and its partial derivatives with respect to the changed coordinate system. In this case, the objective is to find a solution that agrees at the center of the domain, when approached from various directions, with minimum error.

5 Conclusions

In this paper, we have presented a dynamic programming based formulation for obtaining a best piecewise linear approximation for continuous vector fields with the solution constrained by arbitrary boundary conditions. The method attempts obtain a fast but reasonably accurate approximate solution. The method also assumes a discrete space. The parameters to be determined are of the form $\pi_{i}=(A_{i}(t),\,{\mathbf{b}}_{i}(t),\,\theta_{i})$ used to approximate the vector field in the interval $[t_{i},\,t_{i+1}]$, for $0\leq i\leq N-1$. The components $A_{i}(t)$ and ${\mathbf{b}}_{i}(t)$ contain undetermined linear combinations of some fixed basis functions (such as polynomials) defined on the interval $[t_{i},\,t_{i+1}]$. The objective function for minimization of error can be chosen to be the standard ${\mathscr{L}}^{2}$-norm or it can also be chosen from among a large class of error functionas including the uniform metric. The particular aspect of the error functional that allows a formulation of the dynamic programming method is usually a recurrence relation that an associated $k$-step functional satisfies. In the proposed method, the $k$-step functionals satisfy a one-step recurrence relation. It is possible to formulate dynamic programming based methods also for error functionals satisfying more general recurrence relations, involving difficult parameter dependencies. A method for improving an initial approximation by successive iterations is also presented. The proposed method also updates the boundary values. The correction function in the interior of the interval is found by Newton-Raphson method. The instant when to switch from the dynamic programming method to successive iterations method for improvement of the current solution can be determined based on the width of the convergence set for the Newton-Raphson method. The instant when to switch from the dynamic programming method to successive iterations method for improvement of the current solution can be determined based on the width of the convergence set for the Newton-Raphson method, when it is preferred instead of gradient descent method, for this purpose. It is also possible to consider taking a convex combination of formulated updations, for choosing the actual update. If any one of the updation formula always produces a more accurate solution than the other, then the convex combination degenerates into binary exclusive combination.