Regularization of radial solutions of $p$-Laplace equations, and computations using infinite series

We consider radial solutions of equations with the $p$-Laplace operator in $R^n$. We introduce a change of variables, which in effect removes the singularity at $r=0$. While solutions are not of class $C^2$, in general, we show that solutions are $C^2$ functions of $\displaystyle r^{\frac{p}{2(p-1)}}$. Then we express the solution as an infinite series in powers of $\displaystyle r^{\frac{p}{p-1}}$, and give explicit formulas for its coefficients. We implement this algorithm, using Mathematica software. Mathematica's ability to perform the exact computations turns out to be crucial.


Introduction
Recently there has been an enormous interest in equations with the p-Laplace operator in R n (with p > 1, u = u(x), x ∈ R n ) div |∇u| p−2 ∇u + f (u) = 0 , see e.g., a review by P. Drabek [1], and the important paper of B. Franchi et al [2]. Radial solutions of this equation, with the initial data at r = 0, satisfy ϕ(u ′ ) ′ + n − 1 r ϕ(u ′ ) + f (u) = 0, u(0) = α > 0, u ′ (0) = 0 , (1. 1) where ϕ(v) = v|v| p−2 , p > 1. To guess the form of the solution, let us drop the higher order term and consider n − 1 r ϕ(u ′ ) + f (u) = 0, u(0) = α . (1.2) This is a completely different equation, however, in case p = 2, it is easy to check that the form of solutions is the same: in both cases, it is a series ∞ n=0 a n r 2n (with different coefficients), see P. Korman [4] or [5]. It is natural to guess that the form of solutions will be same for (1.1) and (1.2), in case p = 2 too. With that in mind, let us solve (1.2) in case f (u) = e u . Since α > 0, we see from (1.2) that u ′ (r) < 0 for all r. Then ϕ(u ′ ) = −(−u ′ ) p−1 , and we have Integrating, we get We see that u(r) is a function of r p p−1 , and, for r small, we can expand it as a series u(r) = ∞ n=0 b n r n p p−1 , with some coefficients b n . Motivated by this example, we make a change of variables r → z in (1.1), by letting z 2 = r p p−1 . We expect solutions of (1.1) to be of the form ∞ n=0 c n z 2n , which is a real analytic function of z, if this series converges.
The following lemma provides the crucial change of variables.
Then, for p > 2, the change of variables z 2 = r p p−1 transforms (1.1) into Conversely, if the solution of (1.4) is of the form u = v(z 2 ), with v(t) ∈ C 1 (R + ), then the same change of variables transforms (1.4) into (1.1), for any p > 1.
Then (1.1) becomes which simplifies to which lets us compute u ′′ (0) from the equation (1.4) (the existence of u ′′ (0) is proved later). Indeed, assuming that f (α) > 0, we have smooth, provided that f (u) is smooth. It follows that the solution of p-Laplace problem (1.1) has the form u(r p 2(p−1) ), with smooth u(z). We believe that our reduction of the p-Laplace equation (1.2) to the form (1.4) is likely to find other applications.
We express the solution of (1.1) in the form u(r) = ∞ k=0 a k r k p p−1 , and present explicit formulas to compute the coefficients a k . Interestingly, the coefficient a 1 turned out to be special, as it enters in two ways the formula for other a k . Our formulas are easy to implement in Mathematica, and very accurate series approximations can be computed reasonably quickly. We utilize Mathematica's ability to perform the "exact computations", as we explain in Section 3.

Regularity of solutions in case p > 2
It is well known that solutions of p-Laplace equations are not of class C 2 , in general. In fact, rewriting the equation in (1.1) in the form we see that in case p > 2, u ′′ (0) does not exist. We show that in this case the solution of (1.1) is a C 2 function of r p 2(p−1) . We rewrite the equation in (1.1) as Observe that ϕ −1 (t) = −(−t) We recall the following lemma from J.A. Iaia [3].

Representation of solutions using infinite series
We shall consider an auxiliary problem It follows from the last lemma that any series solution of (3.1) must be of the form ∞ n=0 a n z 2n . The same must be true for the problem (1.4), since for z > 0 it agrees with (3.1). Numerically, we shall be computing the partial sums k n=0 a n z 2n , which will provide us with the solution, up to the terms of order O(z 2k+2 ). Write the partial sum in the form whereū(z) = k−1 n=0 a n z 2n . We regardū(z) as already computed, and the question is how to compute a k . Using the constants defined in (1.3), we let Theorem 3.1 Assume that α > 0, f (u) ∈ C ∞ (R), and f (α) > 0. The solution of the problem (1.4) in terms of a series of the form ∞ k=0 a k z 2k is obtained by taking a 0 = α, then and for k ≥ 2, we have (the following limits exist) n=0 a n z 2n is the previously computed approximation, and To check the accuracy of this computation, we denoted by q(z) the left hand side of (4.2) (with u(z) being the above polynomial of tenth degree), and asked Mathematica to expand q(z) into series about z = 0. Mathematica returned: q(z) = O(z 121 10 ). We have performed similar computations, with similar results, for other values of u(0) = α. For larger values of α, e.g., for α = 2, the computations take longer, but no more than several minutes.
We have obtained similar results for all other f (u) and p that we tried (including the case 1 < p < 2). We wish to stress that in all computations, when the solution of (1.4) was computed up to the order z 2n , the defect function q(z) was at least of order O(z 2n+2 ) near z = 0. This heuristic result is consistent with the Theorem 3.1, but does not seem to follow from it.