# How to approximate this bell-shaped integral function using Mathematica?

Consider the following function : $$\tag{1} f(v) = v^{\frac{d}{2}} \int_v^{\infty} u^{a \,-\, \smash{\frac{d}{2}} \,-\, 1} \; e^{-\, a \, u} \; du,$$ where $a$ and $d$ are two positive constants (parameters). Notice the lower limit of the integral : $v$ is a variable with the constaint $v \ge 0$. In Mathematica, the plot of this function (for $0 \le v < \infty$) shows a bell-shaped curve. I'm giving the Mathematica code here to show it :

function[v_, a_, d_] := v^(d/2)NIntegrate[u^(a - d/2 - 1)Exp[- a u], {u, v, Infinity}]

Plot[{
function[v, 1/4, 2],
function[v, 1/2, 3],
function[v, 1, 4]
}, {v, 0, 10},
Frame -> True,
PlotStyle -> {
Directive[Blue, Thick],
Directive[Red, Thick],
Directive[Orange, Thick]
}
]


Now, I would like to define an approximation of this function around the maximal value $f_{\text{max}} = f(v_0)$ (for an unknown $v_0$), since I can't get an explicit expression from the integral. I guess the approximation should look like some Gauss-like function, but I'm not sure.

What is the best approximation of function (1) around its maximal value ? How to find it using Mathematica ?

EDIT : What I'm looking for is an explicit analytical function of all variables ; $v$, $a$, $d$. Not numerical evaluation. Something that may look like this : $$\tag{2} f(v) \approx A(a, d) \, v^{B(a, d)} \, e^{-\, (v \,-\, C(a, d))^2/D(a, d)},$$ or a Taylor expansion around the maximal value of the function : $$\tag{3} f(v) \approx A(a, d) + B(a, d)(v - C(a, d)) + \frac{1}{2} \, D(a, d)(v - C(a, d))^2,$$ where $A(a, d)$, $B(a, d)$, $C(a, d)$ and $D(a, d)$ are analytical functions of $a$ and $d$. I'm especially interested in the function $C(a, d)$ (position of the max value of the $f$ function).

Take note that the derivative of function (1), set to 0, give this relation : $$\tag{4} f_{\text{max}} = \frac{2}{d} \; v_0^a \; e^{-\, a \, v_0},$$ where $v_0 \equiv v_0(a, d)$ is the position of the max value of (1).

• I realized that my question wasn't specifying clearly that I was looking for a function of all variables ; v, a and d. Sorry. Please, see the edit. – Cham Nov 5 '16 at 18:32
• The coefficient $B(a,d)=0$ in $(3)$ up to the Fermat theorem. – user64494 Nov 5 '16 at 18:44
• @user64494, ah ! Yes, of course ! Since it's an approximation around the max value (in the case of the Taylor expansion). But not in the case of the gaussian approximation. – Cham Nov 5 '16 at 18:45

I'm not exactly sure what you are looking for, but here is my best guess. There does appears to be a closed form expression for your function. In particular,

v^(d/2)Integrate[u^(a-d/2-1)Exp[-a u],{u,v,Infinity},Assumptions->Re[v]>0&&Re[a]>0]


Evaluates to

a^(1/2 (-2 a + d)) v^(d/2) Gamma[a - d/2, a v]


Finding the maximum is as simple as differentiating and setting the expression equal to zero. Writing a function

D[a^(1/2 (-2 a+d)) v^(d/2) Gamma[a-d/2,a v],v]
(*-a^(1+1/2 (-2 a+d)) E^(-a v) v^(d/2) (a v)^(-1+a-d/2)
+1/2 a^(1/2 (-2 a+d)) d v^(-1+d/2) Gamma[a-d/2,a v]*)

maxvalue[a_?NumberQ,d_?NumberQ]:=
Re@FindRoot[-a^(1+1/2 (-2 a+d)) E^(-a v) v^(d/2) (a v)^(-1+a-d/2)
+1/2 a^(1/2 (-2 a+d)) d v^(-1+d/2) Gamma[a-d/2,a v]==0,{v,.1}][[1,2]]


Evaluating maxvalue gives the value $v$ for which $f\left(v,a,d\right)$ is maximized. It is also easy to calculate the corresponding maximum value of $f\left(v,a,d\right)$. For example,

maxvalue[1/4, 2]
(*0.321069*)

function[maxvalue[1/4, 2], 1/4, 2]
(*0.694688*)


EDIT

Alternatively, you could solve the problem numerically if a closed form solution to the integral does not exist. In this case I used ND from the NumericalCalculus package.

Needs["NumericalCalculus"]
maxvaluenum[a_?NumberQ,d_?NumberQ]:=FindRoot[ND[function[v,a,d],v,k]==0,{k,.1}][[1,2]]
maxvaluenum[1/4, 2]
(*0.321067*)


This gives roughly the same result as before. I should note that this method is a lot more computationally intensive then the direct method.

EDIT 2

Here is my attempt at what I think you are trying to get at. The general idea is to determine the best fitting polynomial through some range near the neighborhood of the maximum value of $f$. Here we make a new function

approxfun[orderin_,ain_,din_,rangein_]:=
Module[{order=orderin,a=ain,d=din,range=rangein},
maxval=maxvalue[a,d];
data=Table[{v,Re[a^(1/2 (-2 a+d)) v^(d/2) Gamma[a-d/2,a v]]},
{v,maxval-range,maxval+range,range/(10*order)}];
lm=LinearModelFit[data,Evaluate[Table[x^i,{i,1,order}]],x];
Clear[data];
lm];


Using this function creates a FittedModel of a specified order through the given range around the function's maximum value. Example usage:

approxfun[4, 1/4, 2, .1]
(*FittedModel[0.550941 +1.27435 x-4.10958 x^2+5.8875 x^3-3.44731 x^4]*)


You can see that we have created a best fit polynomial function of order 4 with $v \in v_{f_{max}} \pm range$. Plotting with respect to $f$

You can also directly use Series

Plot[Evaluate[Normal@Series[a^(1/2 (-2 a + d)) v^(d/2) Gamma[a - d/2, a v],
{v, maxvalue[a, d], 4}] /. {a -> 1/4, d -> 2}],
{v, Evaluate[maxvalue[a, d] /. {a -> 1/4, d -> 2}] - .1,
Evaluate[maxvalue[a, d] /. {a -> 1/4, d -> 2}] + .1}]


• Is it possible to get an analytical expression for the maxvalue of the function and its argument ? I'm not looking for numerical evaluations. What I'm looking for is an expression like $$f(v) \approx \text{some simple and regular function of v},$$ around the max value. – Cham Nov 5 '16 at 17:16
• @Cham There's no definite answer to that. You could use, e.g. Series around the maximum, or sample the function and use Interpolation, or fit a parabola near the maximum... (The locations of the maximum are easily found with Marchi's solution.) We can help you with implementation of the method you choose, but "what is the best approximation of function (1)" is ambiguous, and most likely beyond the scope of this site. – corey979 Nov 5 '16 at 17:41
• @corey979, I'm looking for a simple analytical approximation (no numerical evaluations). A Taylor serie would do (around the max value), but a gaussian would probably be better. Since the question may be ambiguous, many answers are allowed. – Cham Nov 5 '16 at 17:44
• @corey979 Beat me to the punch while working on updating my answer. I was thinking along the same lines, now seeing your answer, obviously. – Marchi Nov 5 '16 at 17:59

Here is another approach by making use of the PadeApproximant around the same point $v=0.5$ for each $d \in [1,4]$ and each $a \in [1/4,2]$. One does not need to find any extrema of the function under consideration.

function[v_, a_, d_] :=  v^((1/2)*d)*NIntegrate[u^(a - (1/2)*d - 1)*Exp[-a*u], {u, v, Infinity},Method -> "DoubleExponential"];

a=PadeApproximant[function[v, 1/4, 2], {v, 0.5, {3, 3}}]


(0.681942 + 1.88855 (-0.5 + v) + 1.06914 (-0.5 + v)^2 - 0.1788 (-0.5 + v)^3)/(1.00000000000000 + 2.94577 (-0.5 + v) + 2.35946 (-0.5 + v)^2 + 0.433973 (-0.5 + v)^3)

Plot[{a, function[v, 1/4, 2]}, {v, 0.1, 2}]


• I edited the question to make a precision about what I'm looking for. I hope it's clearer. – Cham Nov 5 '16 at 18:27

Taking function from the OP

function[v_, a_, d_] :=
v^(d/2) NIntegrate[u^(a - d/2 - 1) Exp[-a u], {u, v, Infinity}]


and maxvalue from Marchi's answer

maxvalue[a_?NumberQ, d_?NumberQ] :=
Re@FindRoot[-a^(1 + 1/2 (-2 a + d)) E^(-a v) v^(d/2) (a v)^(-1 + a -
d/2) + 1/2 a^(1/2 (-2 a + d)) d v^(-1 + d/2) Gamma[a - d/2,
a v] == 0, {v, .1}][[1, 2]]


for some parameters

{a, d} = {5/7, 2.4};


I find the max value

max = maxvalue[a, d]


0.453786

Then I sample the function near the maximum:

data = Table[{x, function[x, a, d]}, {x, max - 1/10 max, max + 1/10 max, 1/100 max}];

plot1 = ListPlot @ data;


and fit a parabola

nlm = Normal @ NonlinearModelFit[data, A x^2 + B x + c, {A, B, c}, x]


0.261741 + 0.356048 x - 0.391348 x^2

Verification:

Maximize[nlm, x]


{0.342724, {x -> 0.4549}}

plot2 = Plot[nlm, {x, max - 1/10 max, max + 1/10 max}, PlotStyle -> Red];


Visually the approximation is very satisfying:

Show[plot1, plot2]
`

One can try with different functions (maybe a Gaussian curve), different range of the region to approximate in, etc.

• This is still numerical. I strongly suspect (I may be wrong !) that we could find an analytical function with parameters (i.e. symbolic function) : $$f(v) \approx g(v, a, d),$$ where the g function is a gaussian or a polynomial. – Cham Nov 5 '16 at 18:13
• Please, see the edited question. – Cham Nov 5 '16 at 18:27