Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a non-linear integral equation that I'd like to solve with Mathematica:

$$ \int_{0}^{1} \mathrm{d}x \frac{B(x) v}{(B(x) + B(v))^2} = 1$$

Solve[Integrate[B[x] v / (B[x] + B[v])^2, {x, 0, 1}] == 1, B[v], {v}]

This problem comes up in an exotic problem from auction theory where I want to find the optimal bidding function B[v] for v in [0,1]. Is there any way to get a solution to this problem in Mathematica? Closed-form solutions are preferred (duh!) but a plot of a numerical one is also OK.

share|improve this question
2  
Mathematica doesn't have any functionality for solving Integral equations. They must, if possible, be first converted to differential equations. You could numerically solve the equation by putting it in the form B'[v] = f[v] where f[v] is a function of v that has an integral somewhere in it, but that doesn't look possible in this case. –  Searke Mar 11 '13 at 18:44
    
@Searke OK, thanks for the info. –  TemplateRex Mar 11 '13 at 19:06

1 Answer 1

up vote 18 down vote accepted

Let me show how to roll your own numerical solution to a non-linear integral equation using a collocation method. It's fun!

This will involve two approximations. First, we will approximate the function B[x] by its values at n particular points in the range {x, 0, 1}. The integral over x will be replaced by a weighted sum over n, i.e., a quadrature rule. Second, we will only exactly satisfy the integral equation at those n points. It will hopefully be approximately satisfied at other points.

Quadrature rule

Let's borrow one of NIntegrate's quadrature rules:

order = 5;
{abscissae, weights, errweights} = 
 NIntegrate`GaussKronrodRuleData[order, MachinePrecision]

enter image description here

The first part abscissae is the list of n points in the range {0, 1} at which we will approximate the solution. The second part weights is the vector of weights for function values at those points to compute an integral estimate.

n = Length[abscissae]

enter image description here

Approximate integral equation

Using the weights of the quadrature rule we can define the approximate integral equation at a single point:

integralEquation[vals_, v_, fv_] := -1 + 
  weights.Table[fx v/(fx + fv)^2, {fx, vals}]

Here vals is the list of function values at the abscissae, representing the function B[x] in your question, and v and fv are particular values of v and B[v] in your question.

The value of the above function is zero when the integral equation is approximately satisfied at the specified point.

Here is a vector-valued version that evaluates the approximate integral equation at all of the abscissae:

integralEquation[vals_] := 
 MapThread[integralEquation[vals, ##] &, {abscissae, vals}]

Given the vector of function values at the absissae, it gives the vector of residuals indicating whether the integral equation is satisified:

integralEquation[RandomReal[1, n]]

enter image description here

We would like to find vals such that the residuals are all zero! I.e., we want to find a root of integralEquation.

Find root

Let's use the symbols {c[1], ..., c[n]} to represent the solved function values at the abscissae:

vars = Array[c, n]

enter image description here

We will find a root using FindRoot, which likes to have a decent guess for the values of the c[i] to start with. By experimenting a little I learned that a linear solution of the form B[x]==0.2*x is a decent guess:

guess = Thread[{vars, 0.2 abscissae}]

enter image description here

Now we can search for a root of integralEquation:

solution = FindRoot[integralEquation[vars], guess]

enter image description here

This represents an approximate solution to the integral equation.

Solution function

We can plot the approximate solution over the whole range {0, 1} by interpolating the solution we just obtained at the abscissae:

f[x_] = InterpolatingPolynomial[Thread[{abscissae, vars}] /. solution, x];
Plot[f[x], {x, 0, 1}]

enter image description here

The solution is not exactly correct. In fact, for very small values of x, the residual (the amount by which the integral equation is not satisfied) is substantial. Plot the residual for a few values of x:

residual[f_, v_] := 
 Module[{x}, -1 + NIntegrate[f[x] v/(f[x] + f[v])^2, {x, 0, 1}]];
ListLinePlot[Table[{x, residual[f, x]}, {x, 0.01, 1, 0.01}], 
 Frame -> True, PlotRange -> All, AxesStyle -> Dashed]

enter image description here

From here you have some options.

  • You could increase the value of n (by increasing the value of order). You may run into issues with working precision (machine precision may not be enough). You could also use a local interpolation (Interpolation) instead of a global interpolation (InterpolatingPolynomial), which might be more stable for large n.
  • You might be able to solve the tricky bit near x==0 with other methods such as a power series expansion.

Here is a notebook containing the above prototype.

Linear integral equation

For the benefit of people solving related problems, let me just mention that we used FindRoot to search for a root (starting from a plausible guess) because this is a nonlinear integral equation. For a linear integral equation, you can use the same collocation method, but the integralEquations will be linear, so you can be sure of finding a solution simply by using Solve, or reformulate as a matrix problem and use LinearSolve.

share|improve this answer
    
This is a great answer. +1 –  gpap Mar 12 '13 at 18:30
    
Wow... A really awesome answer, never expected such a level of detail! +1 and accepted. I think the x=0 irregularity could be resolved by imposing B[x]==0 (i.e. never bid for things one values at zero). If even that fails, B'[x]>=0 (i.e. bids are monotonically increasing in the valuation) might give additional constraints on interpolating functions. Could you give a hint on how to impose such boundary conditions? –  TemplateRex Mar 12 '13 at 19:52
1  
You could impose B[x]==0 as one of the n collocation points (a special case would be required for integralEquation to avoid appearing to divide by 0). A collocation point at the other endpoint would be fine too. Try LobattoKronrodRule for quadrature rules that include endpoints. –  Andrew Moylan Mar 12 '13 at 20:34
    
I will try that. Thanks again for this wonderful answer. I learnt a lot! –  TemplateRex Mar 12 '13 at 20:54
    
@AndrewMoylan Do you think the strategy suggested here of increasing the point's density near the interval endpoints could be of use here? –  belisarius Mar 12 '13 at 22:40

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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