Quasilinearization

Two numerical solutions of the nonlinear example boundary value problem ${\displaystyle y''=y^{2}}$, ${\displaystyle y(-1)=y(1)=1}$. Solved by a spectral Chebyshev method and quasilinearization. The top curve ${\displaystyle u_{1}}$ used 21 interpolation nodes, and the bottom curve ${\displaystyle u_{2}}$ used 34. Both used 3 iterations.

In mathematics, quasilinearization is a technique which replaces a nonlinear differential equation or operator equation (or system of such equations) with a sequence of linear problems, which are presumed to be easier, and whose solutions approximate the solution of the original nonlinear problem with increasing accuracy. It is a generalization of Newton's method; the word "quasilinearization" is commonly used when the differential equation is a boundary value problem. ^{[1] [2]}

Abstract formulation

Quasilinearization replaces a given nonlinear operator N with a certain linear operator which, being simpler, can be used in an iterative fashion to approximately solve equations containing the original nonlinear operator. This is typically performed when trying to solve an equation such as N(y) = 0 together with certain boundary conditions B for which the equation has a solution y. This solution is sometimes called the "reference solution". For quasilinearization to work, the reference solution needs to exist uniquely (at least locally). The process starts with an initial approximation y₀ that satisfies the boundary conditions and is "sufficiently close" to the reference solution y in a sense to be defined more precisely later. The first step is to take the Fréchet derivative of the nonlinear operator N at that initial approximation, in order to find the linear operator L(y₀) which best approximates N(y)-N(y₀) locally. The nonlinear equation may then be approximated as N(y) = N(y_k) + L(y_k)( y - y_k) + O( y-y_k )², taking k=0. Setting this equation to zero and imposing zero boundary conditions and ignoring higher-order terms gives the linear equation L(y_k)( y - y_k ) = - N(y_k). The solution of this linear equation (with zero boundary conditions) might be called y_k+1. Computation of y_k for k=1, 2, 3,... by solving these linear equations in sequence is analogous to Newton's iteration for a single equation, and requires recomputation of the Fréchet derivative at each y_k. The process can converge quadratically to the reference solution, under the right conditions. Just as with Newton's method for nonlinear algebraic equations, however, difficulties may arise: for instance, the original nonlinear equation may have no solution, or more than one solution, or a multiple solution, in which cases the iteration may converge only very slowly, may not converge at all, or may converge instead to the wrong solution.

The practical test of the meaning of the phrase "sufficiently close" earlier is precisely that the iteration converges to the correct solution. Just as in the case of Newton iteration, there are theorems stating conditions under which one can know ahead of time when the initial approximation is "sufficiently close".

Contrast with discretizing first

One could instead discretize the original nonlinear operator and generate a (typically large) set of nonlinear algebraic equations for the unknowns, and then use Newton's method proper on this system of equations. Generally speaking, the convergence behavior is similar: a similarly good initial approximation will produce similarly good approximate discrete solutions. However, the quasilinearization approach (linearizing the operator equation instead of the discretized equations) seems to be simpler to think about, and has allowed such techniques as adaptive spatial meshes to be used as the iteration proceeds. ^[3]

Example

As an example to illustrate the process of quasilinearization, we can approximately solve the two-point boundary value problem for the nonlinear node $${\displaystyle {\frac {d^{2}}{dx^{2}}}y(x)=y^{2}(x),}$$ where the boundary conditions are ${\displaystyle y(-1)=1}$ and ${\displaystyle y(1)=1}$. The exact solution of the differential equation can be expressed using the Weierstrass elliptic function ℘, like so: ${\displaystyle y(x)=6\wp (x-\alpha |0,\beta )}$ where the vertical bar notation means that the invariants are ${\displaystyle g_{2}=0}$ and ${\displaystyle g_{3}=\beta }$. Finding the values of ${\displaystyle \alpha }$ and ${\displaystyle \beta }$ so that the boundary conditions are satisfied requires solving two simultaneous nonlinear equations for the two unknowns ${\displaystyle \alpha }$ and ${\displaystyle \beta }$, namely ${\displaystyle 6\wp (-1-\alpha |0,\beta )=1}$ and ${\displaystyle 6\wp (1-\alpha |0,\beta )=1}$. This can be done, in an environment where ℘ and its derivatives are available, for instance by Newton's method. ^{[lower-alpha 1]}

Applying the technique of quasilinearization instead, one finds by taking the Fréchet derivative at an unknown approximation ${\displaystyle y_{k}(x)}$ that the linear operator is $${\displaystyle L(\varepsilon )={\frac {d^{2}}{dx^{2}}}\varepsilon (x)-2y_{k}(x)\varepsilon (x).}$$ If the initial approximation is ${\displaystyle y_{0}(x)=1}$ identically on the interval ${\displaystyle -1\leq x\leq 1}$, then the first iteration (at least) can be solved exactly, but is already somewhat complicated. A numerical solution instead, for instance by a Chebyshev spectral method using ${\displaystyle n=21}$ Chebyshev—Lobatto points ${\displaystyle x_{k}=\cos(\pi (n-1-k)/(n-1))}$ for ${\displaystyle k=0,1,\cdots ,n-1}$ gives a solution with residual less than ${\displaystyle 5\cdot 10^{-9}}$ after three iterations; that is, ${\displaystyle y_{3}(x)}$ is the exact solution to ${\textstyle {\frac {d^{2}}{dx^{2}}}y(x)-y^{2}(x)=5\cdot 10^{-9}v(x)}$, where the maximum value of ${\displaystyle |v(x)|}$ is less than 1 on the interval ${\displaystyle -1\leq x\leq 1}$. This approximate solution (call it ${\displaystyle u_{1}}$) agrees with the exact solution ${\displaystyle 6\cdot \wp (x-\alpha |0,\beta )}$ with $${\displaystyle \{\alpha \approx 3.524459420,\beta \approx 0.006691372637\}.}$$

Other values of ${\displaystyle \alpha }$ and ${\displaystyle \beta }$ give other continuous solutions to this nonlinear two-point boundary-value problem for ODE, such as $${\displaystyle \{\alpha \approx 2.55347391110,\beta \approx -1.24923895273\}.}$$ The solution corresponding to these values plotted in the figure is called ${\displaystyle u_{2}}$. Yet other values of the parameters can give discontinuous solutions because ℘ has a double pole at zero and so ${\displaystyle y(x)}$ has a double pole at ${\displaystyle x=\alpha }$. Finding other continuous solutions by quasilinearization requires different initial approximations to the ones used here. The initial approximation ${\displaystyle y_{0}=5x^{2}-4}$ approximates the exact solution ${\displaystyle u_{2}}$ and can be used to generate a sequence of approximations converging to ${\displaystyle u_{2}}$. Both approximations are plotted in the accompanying figure.

Notes

1. For more information about elliptic functions, see Lawden (1989). ^[4]

See also

• Describing function

References

1. Ascher, Uri M.; Mattheij, Robert M.; Russell, Robert D. (1995). Numerical solution of boundary value problems for ordinary differential equations. SIAM.
2. Sylvester, R.J.; Meyer, F. (June 1965). "Two point boundary problems by quasilinearization". Journal of the Society for Industrial and Applied Mathematics. 13 (2): 586–602. doi:10.1137/0113038. JSTOR   2946451 . Retrieved 31 March 2022.
3. Bornemann, Folkmar (1991). An Adaptive Multilevel Approach to Parabolic Equations in Two Space Dimensions. ZIB. Retrieved 6 March 2022.
4. Lawden, Derek F. (1989). Elliptic Functions and Applications. New York: Springer-Verlag. ISBN   0-387-96965-9.

Further reading

• https://encyclopediaofmath.org/wiki/Quasi-linearization