I'm trying to approximate a decaying at infinity 2D function f[x_,y_] with a rational function of 2 variables. For a 1D function like that PadeApproximant would be a natural choice. Unfortunately, documentation seems to indicate a 1D PadeApproximant only. Is is possible to approximate a 2D function with a rational function?


There are known approximations of integral over an interval of 1D Gaussian distribution. It's elementary to extend that to an integral over a rectangle of a 2D Gaussian distribution. I am trying to approximate an integral over a triangle of a 2D Gaussian. This boils down to approximating Integral[Erf[ax+b]*E^(-x^2),x]. Asking Mathematica directly that yields nothing. So I'm looking for the ways to make Mathematica fit that expression with elementary functions.

  • 1
    $\begingroup$ Good idea to supply the function that you are attempting to approximate. What have you tried so far? $\endgroup$ – Jack LaVigne Aug 29 '16 at 19:13
  • $\begingroup$ Remember that Pade Approximations are most accurate near about a particular point and can be singular far from that point. It may be better to use FindFit[] or similar functions to obtain a rational approximation that is reasonably accurate over a larger domain, although less accurate than a Pade Approximation near a particular point. $\endgroup$ – bbgodfrey Aug 29 '16 at 22:15
  • $\begingroup$ @bbgodfrey: What would be the symbolic versions of FinFit? I'm aware of Series, which would give me Taylor series (in any dimension), and that is asymptotically wrong for me because my function decays at infinity, and of PadeApproximant, which seems to work only in 1D. Does Mathematica has other symbolic approximations, either in rational functions or in elementary functions? Fourier series won't work for me either because finite number of terms would sum up to a periodic function, and mine decays at infinity. $\endgroup$ – Michael Aug 30 '16 at 0:01
  • 1
    $\begingroup$ It seemed fairly simple to give a Chebyshev series approximation over a bounded interval, but you seemed uninterested... $\endgroup$ – Michael E2 Sep 1 '16 at 12:11
  • 1
    $\begingroup$ The EDIT is close to a duplicate of mathematica.stackexchange.com/questions/121963/… $\endgroup$ – Michael E2 Sep 1 '16 at 12:14

Perhaps surprisingly, the Pade Approximant for the quantity,

Integrate[Erf[x + b]*E^(-x^2), {x, 0, c}];

can be obtained quite easily.

pol = PadeApproximant[Integrate[Erf[x + b]*E^(-x^2), {x, 0, c}], {c, 0, 6}] // Simplify;

Admittedly, this is a 1D expansion, not a 2D expansion as requested in the question. However, the resulting expression is analytical, so there seems little point in expanding it in b. Note that the LeafCount is large, 5911.

An accurate numerical expression is needed for comparison, and the method I have chosen is

sol = ParametricNDSolveValue[{y'[x] == Erf[x + b]*E^(-x^2), y[0] == 0},
    y, {x, 0, c}, {b, c}];

The comparison,

Plot3D[{sol[b, c][c], pol}, {b, -2, 2}, {c, 0, 4}, 
    AxesLabel -> {b, c, s}, ViewPoint -> {-3, 8, 3}]

enter image description here

indicates superb agreement for c less than about 2, but rather less good agreement for larger values, especially when b is negative. Note, however, that the integral is nearly constant for c greater than about 2.5, and that value could be used there in place of the PadeApproximant. (The asymptotic value of the integral would be computed numerically.)

It should also be noted that the integral can be integrated by parts to yield,

Erf[c] Erf[b + c] Sqrt[Pi]/4 + 
    Integrate[(Erf[x + b]*Exp[-x^2] - Erf[x]*Exp[-(x + b)^2])/2, {x, 0, c}]

A particular advantage of this representation is that the integral here vanishes for b == 0, leaving

Erf[c]^2 Sqrt[Pi]/4.

The expansion of the integral in this case is

pol2 = PadeApproximant[Integrate[(Erf[x + b]*Exp[-x^2] - Erf[x]*Exp[-(x + b)^2])/2, 
    {x, 0, c}], {c, 0, 6}] // Simplify;

Its LeafCount is modestly smaller, 5033. Comparison with numerical results shows that this representation is better for some values of b (especially those near zero) but worse for some other values.

Plot3D[{sol[b, c][c], Erf[c] Erf[b + c] Sqrt[Pi]/4 + pol2}, {b, -2, 2}, {c, 0, 4},
    AxesLabel -> {b, c, s}, ViewPoint -> {-3, 8, 3}]

enter image description here

Other Mathematica functions exist that could be used to obtain other rational approximations to the integral, such as RationalInterpolation, EconomizedRationalApproximation, and even FindFit. An interesting approach is provided in the article, Generalized Padé Approximation. I have not experimented with these other approaches.


Expanding Erf[x + b] as a Taylor series and inserting it into the integral also works well for positive b, but definitely not for negative b.

Series[Erf[x + b], {x, b, 12}] // Normal;
tol = Integrate[s*E^(-x^2), {x, 0, c}, Assumptions -> c ∈ Reals // Simplify;
Plot3D[{sol[b, c][c], tol}, {b, -2, 2}, {c, 0, 4}, 
    AxesLabel -> {b, c, "s"}, ViewPoint -> {-3, 8, 3}]

enter image description here

LeafCount is 898.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.