An Efficient Approach for Solving a Class of Nonlinear 2d Parabolic Pdes

We consider a class of nonlinear 2D parabolic equations that allow for an efficient application of an operator splitting technique and a suitable linearization of the discretized problem. We apply our scheme to study the finite extinction phenomenon for the porous-medium equation with strong absorption. A comparison of accuracy and computational efficiency of the resulting algorithms for several test problems is presented. 1. Introduction. We study a class of nonlinear 2D parabolic PDEs where the non-linearity is a power of the solution. We apply a linearization and an operator splitting technique. We use our algorithm for computing to high accuracy the extinction time for the porous-medium equation with strong absorption. We use a finite-difference method and an implicit scheme of the Crank-Nicolson type in view of its superior stability properties. Then we are led to systems of nonlinear algebraic equations. These can be solved using Newton's type methods, which is costly. Instead, we choose to use a linearization approach that eliminates the need for iterations at each temporal step. Solving the arising large linear system of algebraic equations straightforwardly using direct methods would be extremely inefficient. Thus, the linearization by itself does not solve the main difficulty, the high computational cost, in dealing efficiently with 2D problems. To remedy this, to the linearized equations we apply the operator splitting technique that allows for computationally more efficient solution processes. The nonlinearity in the considered PDEs has the form of powers of the solution function. This is a quite common situation in many applications, especially, in the porous-medium type equations, see [1, 4], for which various numerical schemes have been introduced, see, for example, [4, 6, 7]. We apply our scheme in the porous-medium equation , and we use it to study the finite-extinction phenomenon for the porous-medium equation with strong absorption. Throughout the paper, we use the following finite-difference operator notation. A continuous function u = u(x, y, t) is discretized on an equidistant spatio-temporal grid (x i ,y j ,t n) = (ih, jh, nk), 1 ≤ i, j ≤ M, 0 ≤ n ≤ N, where h = ∆x = ∆y and k = ∆t are step sizes in spatial and temporal directions and n = 0 is for an initial function. The discrete function at the point (x i ,y j ,t n) is denoted by u i,j,n. We define U n as a discrete column vector of …


Introduction.
We study a class of nonlinear 2D parabolic PDEs where the nonlinearity is a power of the solution.We apply a linearization and an operator splitting technique.We use our algorithm for computing to high accuracy the extinction time for the porous-medium equation with strong absorption.
We use a finite-difference method and an implicit scheme of the Crank-Nicolson type in view of its superior stability properties.Then we are led to systems of nonlinear algebraic equations.These can be solved using Newton's type methods, which is costly.Instead, we choose to use a linearization approach that eliminates the need for iterations at each temporal step.Solving the arising large linear system of algebraic equations straightforwardly using direct methods would be extremely inefficient.Thus, the linearization by itself does not solve the main difficulty, the high computational cost, in dealing efficiently with 2D problems.To remedy this, to the linearized equations we apply the operator splitting technique that allows for computationally more efficient solution processes.
The nonlinearity in the considered PDEs has the form of powers of the solution function.This is a quite common situation in many applications, especially, in the porousmedium type equations, see [1,4], for which various numerical schemes have been introduced, see, for example, [4,6,7].We apply our scheme in the porous-medium equation, and we use it to study the finite-extinction phenomenon for the porous-medium equation with strong absorption.
Throughout the paper, we use the following finite-difference operator notation.A continuous function u = u(x, y, t) is discretized on an equidistant spatio-temporal grid (x i ,y j ,t n ) = (ih, jh, nk), 1 ≤ i, j ≤ M, 0 ≤ n ≤ N, where h = ∆x = ∆y and k = ∆t are step sizes in spatial and temporal directions and n = 0 is for an initial function.
The discrete function at the point (x i ,y j ,t n ) is denoted by u i,j,n .We define U n as a discrete column vector of order M 2 at the time level n: U n = u 1,1,n ,u 2,1,n ,...,u M,M,n T . (1.1) The forward difference operators of each direction are defined as D +y : D +y u j = u j+1 − u j h . (1.2) To simplify the notation, indices over variables that do not change under the particular operator (i.e., temporarily frozen) are suppressed.The backward difference operator is defined analogously as These one-sided operators are first-order accurate approximations of the derivatives, that is, (D +x −D)u i = ᏻ(h), and so forth.Additionally, we use the second-order operator which is second-order accurate, that is, Note that the result of these operators is a column vector or a matrix depending on contexts when the operators are applied to column vectors.This paper is arranged as follows.In Section 2, we present our model problem.Section 3 is devoted to the analysis of computational efficiency of the different algorithms for porous-medium equation.Numerical experiments for several problems with various conditions are reported in Section 4, while Section 5 contains the concluding remarks.

Model problem: porous-medium equation.
We consider the following initialboundary value problem: ) where ∆ denotes the Laplace operator, m ≥ 2, p > 0, c is a constant, the bounded spatial region , and I = (0, T), T < ∞, is a time interval.We let g and u 0 be given smooth data.We will use problem (2.1) as a model apt to describe our method.
If m ≥ 2 and c = 0, (2.1a) is a well-known porous-medium equation which governs the density of an ideal gas [1].If m ≥ 2, c > 0, and 0 < p < 1, (2.1a) is called the porous-medium equation with strong absorption.It appears in various applications, most frequently to describe a phenomenon of thermal propagation in an absorptive medium, where the solution u stands for temperature.In other applications, u is a concentration and the process is described as diffusion with absorption.The solution u is required to be nonnegative, and thus we consider nonnegative and bounded initial condition u 0 (x, y).Furthermore, in the case of zero boundary conditions, that is, g(x) = 0, it is known (see [2]) that the solution u vanishes identically as time t goes on and extinction in finite time arises, that is, there exists a time T * > 0 such that u(•,t) is not identically zero for 0 < t < T * but u(•, T * ) ≡ 0. Here, T * is called an extinction time of a solution u and the behavior is known as finite-extinction phenomenon.

Linearization approach. Consider the model porous-medium equation (without absorption)
Applying the standard Crank-Nicolson scheme (which is ᏻ(k 2 + h 2 )-accurate), one obtains where the power m to discrete vector U is done element by element.Rearranging the U n and U n+1 terms gives This yields a system of nonlinear algebraic equations that can be solved using Newton's iterative method at each temporal step.At (n + 1)th temporal step, this system of nonlinear equations is where , and U n is the computed solution vector at the previous time level, a known vector.Then Newton's iterative method is given by J U (l) W (l+1) = −F U (l) , l≥ 0, (3.5)where (l) is an iterative index, W (l+1) = U (l+1) −U (l) , J(U (l) ) is the M 2 ×M 2 Jacobi matrix given by and I is the M 2 × M 2 identity matrix.
The feature of the Crank-Nicolson scheme is its unconditional stability, see the appendix.For all values of λ = k/h 2 , in practice, however, one uses moderate values of λ to avoid slowly decaying oscillations.The main problem one needs to deal with in solving (3.3) is its computational cost.Although one might improve the efficiency of the procedure by exploiting sparsity (see [10]), the Jacobian matrices involved here typically are very large, thus the computational cost is large.

Standard linearization method.
Here, we follow the idea of Richtmyer and Mortan [8] who applied it to 1D problems u t = (u m ) xx (some call it Richtmyer's linearization method, see [9]).Recently, this idea was also applied to Korteweg-de Vries equation, see [5].
For m ≥ 2, or where ω n+1 ≈ U n+1 − U n .Substituting (3.8) into (3.3) and rearranging, we obtain (3.9) Equation (3.9) needs to be solved at each time level n = 1, 2,...,N.In matrix form, it can be written as a linear system in the unknown ω n+1 : where A is M 2 × M 2 block tridiagonal matrix: each of the blocks T j 's and D j 's are M × M tridiagonal and diagonal matrices, respectively.They are not constant diagonal (Toeplitz) but depend on the solution u at the previous time level n: where we can separate the boundary of Ω, ∂Ω, into two parts as follows: (3.15)So we can write U m l | ∂Ω , for l = n, n + 1, the contribution from the boundary: where indices 0 and M + 1 denote the boundary values.The matrix B is an M 2 × M 2 block tridiagonal matrix, I is the M × M identity matrix, and T is an M × M tridiagonal matrix: It should be noted that the linearized equation (3.9) is identical to those obtained by using only one iteration of Newton's method with the initial guess U (0) n+1 = U n , the computed solution vector at the previous time level.This can be verified simply by inspection.The cited papers failed to make this observation.
Again, the straightforward attempt to solve system (3.10) is unacceptably costly.The half bandwidth of A is M, and thus applying a band matrix solver would cost ᏻ(M 2 ) flops per unknown, which is far from optimal.

Operator splitting method.
In order to solve system (3.10),we propose an operator splitting technique.
We add a benign (because it does not affect the ᏻ(k 2 + h 2 )-accuracy of the Crank-Nicolson scheme) term with a truncation error ᏻ(k 2 ): to the only left-hand side of (3.9).As a result we can factorize the left-hand side of (3.9) as follows: n with the suitable boundary contribution.By splitting (3.19) into two simpler systems, we have where the intermediate boundary values can be given from (3.21) by where where for 1 ≤ j ≤ M. For (3.21) the indices in T j are transposed: u m−1 j,i,n instead of u m−1 i,j,n , for 1 ≤ i, j ≤ M. Since the entries depend on the solution u, the matrices in (3.20) and (3.21) are identical only if the solution is symmetric in x and y.
Each tridiagonal block can be solved independently of other blocks at the cost proportional to its size, that is, ᏻ(M).Since there are 2M blocks to be solved, the total cost is ᏻ(M 2 ) or ᏻ(1) flops per unknown, which is optimal.Now we study the finite-extinction phenomenon for (2.1a) with m ≥ 2, 0 < p < 1, and c > 0. Applying the operator splitting technique, we obtain a system of two equations which is similar to the system in (3.20) and (3.21) except b in (3.20 Here, however, in the presence of the additional nonlinear zero-order term cu p , (3.20) becomes nonlinear (because of the term U p n+1 in b).To solve it, we use Newton's method, again with the Jacobian matrix that is block diagonal, and thus each tridiagonal block can be solved independently.As a consequence, the computational cost can significantly be reduced as in the case of the porous-medium equation (without absorption).

Numerical experiments.
In this section, we consider two different numerical experiments.First, we investigate the efficiency of our linearization with splitting method for solving 2D problems.We compare three different methods: two iterations of Newton's method (denoted by N in tables and figures), and the linearization method, first in its standard form (denoted by L), then combined with the operator splitting (denoted by S), in terms of the computational efficiency (cost).We also comment on the accuracy of the schemes.Second, we study the finite-extinction phenomenon for the porous-medium equation with strong absorption employing our scheme.We present the experimental results of the numerical extinction time values for various spatial and temporal step sizes.All our programs, written in Matlab, are implemented on a PC with Pentium III processor at 933 MHz.

Efficiency of the splitting method.
We consider problems defined in a unit square spatial domain divided equidistantly into M grid points in the x-and y-directions.Thus, the spatial step size of the uniform grid becomes h = 1/(M + 1).To study the computational efficiency (cost) of the schemes (per temporal step), we compare the CPU time and the number of Mflops for M = 20, 40, and 80 with a fixed value of λ = k/h 2 = 1, where k is the temporal step size.We report on the cost reduction factor as well as the power, p, of the cost's growth model: cost = ᏻ(M p ) (computed as the slope of the least-squares linear polynomial of logarithmically scaled points).The total cost is then increasing linearly with the number of temporal steps, that is, proportional to the value of 1/λ.It should be noted that the Crank-Nicolson scheme allows for the use of higher values of λ than of the order of 1, see Table 4.1.
The normalized L 2 -norm of the error is defined as where u i,j,n is the numerical solution and u(ih, jh, T ) is the analytical solution at the time t = nk = T , and h and k are the spatial and temporal step sizes, respectively.We investigate two numerical examples which are different initial-boundary value problems of the same porous-medium equation u t = ∆(u 5 ).
Example 4.1.We choose initial and boundary conditions corresponding to the exact solution of u t = ∆(u 5 ) which is defined as (see [10]) for 0 ≤ x, y ≤ 1, and 0 ≤ t ≤ 1.
Example 4.2.Similarly, we consider the initial-boundary value problem corresponding to the exact solution for 0 ≤ x, y ≤ 1, and 0 ≤ t ≤ 0.5.
The exact solution of Example 4.2 (given analytically by Barenblatt, and thus called the Barenblatt solution) is a radially symmetric self-solution u ≥ 0, defined on {(x,t)|x ∈ R N , 0 < t < T } (see [1,3]).It is of the general form The computational cost per temporal step of solving the linear system represented by system (3.10) is the critical component of the total computational expenses.A brute force band Gaussian elimination of the M 2 × M 2 system with the half-band width M would cost ᏻ(M 4 ) flops, a prohibitively high expense.
The pentadiagonal (or block-tridiagonal) matrices (3.11) have the same structure as the discrete Laplacian, although they are not constant diagonal (Toeplitz) and thus do not allow for the use of FFT-based fast solvers.One possibility of reducing the cost is to apply nearly optimal reordering before the elimination.Such a tool is provided in Matlab (as a default option) using the backslash (\) operation in the sparse mode; it employs the minimum degree ordering algorithm.Because of this, the flop count for solving (3.10) decreases to about O(M 3.6 ), see Figure 4.1.At the same time, for the splitting method, see Figure 4.1, the computational cost for solving (3.20) and (3.21) is of about ᏻ(M 2 ) flops, or O(1) per unknown, which is optimal.Thus, our analysis from Section 3.2 is confirmed by the experiments.
One should point out that the CPU growth factors are significantly closer, of about ᏻ(M 3 ) and ᏻ(M 2.5 ), respectively.This measure of efficiency is much more computerand language-dependent.In this implementation, the difference in CPU time between the two methods is less pronounced than the flop count.Nevertheless, for M = 80, the actual cost (in CPU) reduction ratio from Newton's method to our operator splitting method is about 8 times (see Table 4.2).The computational cost for both Examples 4.1 and 4.2 is almost identical, therefore we provide the data only for one of them.
We also provide the comparison of the accuracy of the considered methods.One should remark that the accuracy is heavily problem-dependent as the discretization error depends not only on the numerical scheme and the parameters of discretization but also on the smoothness of the solution.Because of this, we provide the two examples.The computed solution in Example 4.1,   Table 4.1.We also verified the influence of larger values of λ on accuracy, see Table 4.1.
Only for the coarse grids, the slowly decaying oscillations affect the error.
To bring yet another perspective, we plot the comparison of the accuracy versus efficiency of the considered algorithms, see Figure 4.3.It clearly shows the superiority of the splitting method.As a consequence, the results presented in the next section were obtained only by the latter.

Finite-extinction phenomenon.
To study the finite-extinction phenomenon numerically, we consider the initial-boundary value problem (2.1) with m = 2, p = 0.5, c = 5, and zero boundary condition (g = 0).The presented results were obtained by the splitting method.Example 4.3.We choose the initial function u 0 (x, y) given as (see [7]) where We denote M to be the number of spatial steps in the x-and y-directions, that is, (M −1) is the number of grid points in each direction.Thus the step size of the uniform grid is h = 2.4/M.
The second equation of our system, (3.21), is linear.We solve it just as described in Section 4.1.On the other hand, the first equation, (3.20), becomes nonlinear.To solve it we use Newton's iterative method: at the (n+1)th temporal step, the system of (M −1) decoupled nonlinear equations in the x-direction is where for i = 1, 2,...,M − 1, with the known vector of U n and element-by-element operations.
The Jacobian matrix J(U) is where I is the (M − 1) × (M − 1) identity matrix.The Newton iterations are performed until the stopping criterion n+1 < τ, l = 0, 1,..., is satisfied.Here, the initial guess is U (0) n+1 = U n for n = 0, 1,..., and τ is the given tolerance.To prevent the singularity of the Jacobian matrix, we use regularization, that is, J new = J + I, using the machine precision 2.2e−16.Moreover, we force its solution to remain nonnegative at each temporal step.
On the average, the algorithm uses only about two Newton's iterations to solve the nonlinear equation (3.20).As a consequence, the computational cost per temporal step is only about twice that of the cost discussed in Section 4.1 (without the absorption term), see In Figure 4.5 we plot the time evolution of the maximum height, H, of the computed solution, u, of (4.5) in standard and log-log scales.Note that the log of the solution is decreasing rapidly in the vicinity of the extinction time, T * .Additionally, in Figure 4.6 we plot traces of the computed solution, u, on the xy-plane at four different points in its time evolution.
The goal here is to accurately compute the extinction time, T * , the time at which the solution becomes identically zero, u(x, y, T * ) = 0. Thus, one must accurately compute both the solution, u, and the time, T * .The latter task imposes additional restriction: the temporal step k must be smaller or equal to the required precision in determining the value of T * .This point is well illustrated by Table 4.4.
It is clear that the number of significant digits in the computed value of T * is limited by the − log k.It does not mean though that it is sufficient to take small values of k (the temporal step size) as of course the spatial discretization error also plays an important role, see Table 4. 5.
In these experiments with k = 1e − 5, we limited the refining of the spatial grid because of large computational costs.To remedy this, we designed a modified algorithm.From the point of view of accuracy, one could conduct the experiments with much       For a fixed spatial step size, h, the total computational cost is proportional to the number of temporal steps, T * /k for the original algorithm and T o /k o + (T * − T o )/k for the modified one.For example, for h = 0.04, the computed value of T o turned out to be T o = 0.241.In this case, the reduction of the number of temporal steps was from 26 343 to 2 484, or more than tenfold.We have not made an attempt to optimize the speedup.
Although we do not know the exact value of the extinction time, we can obtain a very rough estimate of the power, p, of the asymptotic error model using the Richardson extrapolation method on the values in Table 4.6.We obtain 2 p = (0.26360 − 0.26343)/ (0.26343 − 0.26337) = 17/6.We can thus conclude that the accuracy of this computation is of about ᏻ(h 1.5 ), which is consistent with the experiments in Section 4.1.

Concluding remarks.
For a class of nonlinear parabolic PDEs in 2D rectangular domains, we have constructed an operator splitting algorithm of optimal efficiency.We have verified experimentally for the porous-medium equation that the computational cost of our scheme is O(1) flops per unknown per temporal step while the accuracy remains the same as for two Newton's iterations.
In the presence of an additional zero-order term (strong absorption term), the asymptotic efficiency remains unchanged, O(1), with the leading constant only twice larger.We have modified the algorithm to adaptively change the temporal step size.This allows computing the extinction times extremely accurately and with significant computational savings.

Call for Papers
As a multidisciplinary field, financial engineering is becoming increasingly important in today's economic and financial world, especially in areas such as portfolio management, asset valuation and prediction, fraud detection, and credit risk management.For example, in a credit risk context, the recently approved Basel II guidelines advise financial institutions to build comprehensible credit risk models in order to optimize their capital allocation policy.Computational methods are being intensively studied and applied to improve the quality of the financial decisions that need to be made.Until now, computational methods and models are central to the analysis of economic and financial decisions.
However, more and more researchers have found that the financial environment is not ruled by mathematical distributions or statistical models.In such situations, some attempts have also been made to develop financial engineering models using intelligent computing approaches.For example, an artificial neural network (ANN) is a nonparametric estimation technique which does not make any distributional assumptions regarding the underlying asset.Instead, ANN approach develops a model using sets of unknown parameters and lets the optimization routine seek the best fitting parameters to obtain the desired results.The main aim of this special issue is not to merely illustrate the superior performance of a new intelligent computational method, but also to demonstrate how it can be used effectively in a financial engineering environment to improve and facilitate financial decision making.In this sense, the submissions should especially address how the results of estimated computational models (e.g., ANN, support vector machines, evolutionary algorithm, and fuzzy models) can be used to develop intelligent, easy-to-use, and/or comprehensible computational systems (e.g., decision support systems, agent-based system, and web-based systems) This special issue will include (but not be limited to) the following topics: • Computational methods: artificial intelligence, neural networks, evolutionary algorithms, fuzzy inference, hybrid learning, ensemble learning, cooperative learning, multiagent learning .13) for 1 ≤ j ≤ M and λ = k/h 2 .The M 2 vector b with the contribution from the boundary in (3.10) has the form b

Figure 4 . 1 .
Figure 4.1.Comparison of the efficiency (measured in Mflops and CPU time (ms)) for the Newton (N) method (two iterations), standard linearization (L) method, and our splitting (S) method as a function of the number of grid points, M, for Example 4.1.Here, p is the power of the cost's growth model: cost = ᏻ(M p ).

Figure 4 . 2 .
Figure 4.2.Comparison of the accuracy for the Newton (N) method (two iterations), standard linearization (L) method, and our splitting (S) method as a function of the number of grid points, M, for λ = 1.Here, p is the power of the error's growth model: error = ᏻ(M p ); (a) shows the computed solution in Example 4.1 while (b) shows that in Example 4.2.

Figure 4 . 4 and
Table 4.3.Here, we denote by P the porous-medium equation, Example 4.2, and by A the porous-medium equation with absorption, Example 4.3.Note that in order to make the comparison, we modify Example 4.2 and run it with the same number of the spatial steps, M = 20, 40, 80, and equal exponents, m = 2.

A, p = 2 Figure 4 . 4 .
Figure 4.4.Comparison between the porous-medium equation with strong absorption (A) and the porous-medium equation (P) in terms of (a) the number of Kflops and (b) CPU time (bottom), as a function of the number of spatial steps, M, for λ = 1.Here, p is the power of the cost's growth model: error = ᏻ(M p ).

Figure 4 . 5 .
Figure 4.5.Time evolution of the maximum height, H, of the numerical solution, u, with the spatial step size h = 0.04 and temporal step size k = 0.0001.Plot (b) is the log-log plot of (a).

Figure 4 . 6 .
Figure 4.6.Finite-extinction phenomenon of the numerical solution for Example 4.3 with the spatial step size h = 0.04 and temporal step size k = 0.0001.Here, H is the maximum height of the numerical solution at the time t.

Table 4 .
1. Comparison of the accuracy.

Table 4 .
2. Comparison of the efficiency (average cost per temporal step) for Example 4.1.
the power of the discretization error model: error = ᏻ(M p ) with the theoretical value for smooth enough functions of p = −2.The experiments indicate that the accuracy of our method is almost the same as the one after two iterations of Newton's method, see

Table 4 .
3. Comparison of computational costs (per temporal step) for modified Examples 4.2 and 4.3.(seeTable4.1)thanthose in Table4.5 if it was not for the resolution requirement that k be sufficiently small at the extinction time T * .The time evolution (in the log-log scale) of the maximum height H of the computed solution u, see Figure4.5,

Table 4 .
5. T * for fixed k = 1e-5..The modified algorithm employs a much larger time step, k o for t ≤ T o , and a small k, corresponding to the required resolution, afterwards.In the experiments reported below, see Table4.6, we have chosen an ad hoc value k o = 100k, and the program was adaptively changing the temporal step size at a point when the slope, p, of the maximum height of u becomes p ≤ p o .Again, we have, ad hoc, chosen p o = −20.The reported λ in Table4.6 is that of λ o = k o /h 2 , that is, that for t ≤ T o . *

•
Application fields: asset valuation and prediction, asset allocation and portfolio selection, bankruptcy prediction, fraud detection, credit risk management • Implementation aspects: decision support systems, expert systems, information systems, intelligent agents, web service, monitoring, deployment, implementation