NDSolve to automate shooting method (or others) for free boundary value problem?

Question

I am trying to solve the following system of differential equations using NDSolve:

$$ \begin{align*} f'(s) &= \frac{f(s)}{g(s)-s}\\ g'(s) &= \frac{g(s)}{f(s)-s}\\ f(0) &= 0\\ g(0) &= 0\\ f(\overline{s}) &= 1\\ g(\overline{s}) &= 2 \end{align*} $$

The solution is a Bayes-Nash equilibrium bidding strategy in a first-price auction, as described in a paper by Hubbard and Paarsch.

The right boundary condition is not necessarily known a priori. And there's a singularity at the left boundary, and the system is overdetermined, but solutions do exist.

Actually, this represents a special case with an analytic solution. $\overline{s}=2/3$ and then $f(s) = \frac{2s}{1+\frac{3}{4}s^2}$, $g(s) = \frac{2s}{1-\frac{3}{4}s^2}$. But I want to try to solve it numerically.

One technique Hubbard and Parsch advocate is a shooting technique. Don't specify the left boundary; take a guess for $\overline{s}$ and specify the right boundary, and solve the reverse initial value problem. Then decrease $\overline{s}$ if the solutions diverge, and increase it if they are too far from $f[0]=0,g[0]=0$.

I can do this trial and error by hand and it works. It would be easy to wrap up NDSolve in a root-finding algorithm to do this.

maxBid = 2/3;
soln = NDSolve[{
   {f'[s] == f[s]/(g[s] - s),
    g'[s] == g[s]/(f[s] - s)},
   {f[maxBid] == 1, g[maxBid] == 2}},
  {f, g}, {s, 0, maxBid}] (* gives nearly correct solution *)

However, I'm wondering if there is a way to get NDSolve to handle all this automatically for me. Can the built-in "Shooting" method be used for this? Do any differential equation wizards know of a better way to attack this problem?

Answer 1

Since we should compute maxBid as well, we can consider this problem as optimization problem, and not as BVP. Therefore we start from a functional to be optimized. For this we use colocation method with Haar wavelets (it is my lovely method because it always works) to project equations as follows

A = 0; B = 1; J = 3; M = 2^J; dx = (B - A)/(2*M); 
h1[x_] := Piecewise[{{1, A <= x <= B}, {0, True}}]; 
p1[x_, n_] := (1/n!)*(x - A)^n; 
h[x_, k_, m_] := 
 Piecewise[{{1, 
    Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]}, {-1, 
    Inequality[(1 + 2*k)/(2*m), LessEqual, x, Less, (1 + k)/m]}}, 0]
p[x_, k_, m_, n_] := 
 Piecewise[{{0, x < k/m}, {(-(k/m) + x)^n/n!, 
    Inequality[k/m, LessEqual, x, 
     Less, (1 + 2*k)/(2*m)]}, {((-(k/m) + x)^n - 
       2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, 
         (1 + 2*k)/(2*m) <= 
     x <= (1 + k)/
      m}, {((-(k/m) + x)^n + (-((1 + k)/m) + x)^n - 
       2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, x > (1 + k)/m}}, 0]
xl = Table[A + l*dx, {l, 0, 2*M}]; xcol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2*M + 1}]; 
f1[x_] := Sum[
    af[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   a0*h1[x]; 
f0[x_] := Sum[
    af[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   a0*p1[x, 1] + f10; 
g1[x_] := Sum[
    ag[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   c0*h1[x]; 
g0[x_] := Sum[
    ag[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   c0*p1[x, 1] + g10; 
bc0 = {f0[0] == 0, g0[0] == 0}; 
bc1 = {f0[s1] - 1 == 0, g0[s1] - 2 == 0}; 
var = Flatten[
   Table[{af[i, j], ag[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]; 
cons = Join[bc0, bc1, {0 < s1 < 1}]; 
eqn = {f1[s] - f0[s]/(g0[s] - s), g1[s] - g0[s]/(f0[s] - s)}; eq = 
 Flatten[Table[eqn, {s, xcol}]]; 
varX = Join[{s1, a0, c0, f10, g10}, var];

Here s1 is unknown parameter to be computed. Next step, we minimize equations with constraints

solM1 = NMinimize[{eq.eq, cons}, varX]

Out[]= {8.8844*10^-19, {s1 -> 0.666052, a0 -> 1.14242, 
  c0 -> 8.02142, f10 -> 0, g10 -> 0, af[0, 0] -> 0.545296, 
  ag[0, 0] -> -5.56734, af[0, 1] -> 0.227924, ag[0, 1] -> -0.361955, 
  af[1, 1] -> 0.251853, ag[1, 1] -> -8.14646, af[0, 2] -> 0.0668996, 
  ag[0, 2] -> -0.0743614, af[1, 2] -> 0.149766, ag[1, 2] -> -0.31305, 
  af[2, 2] -> 0.146919, ag[2, 2] -> -1.15702, af[3, 2] -> 0.103255, 
  ag[3, 2] -> -9.6107, af[0, 3] -> 0.0174592, ag[0, 3] -> -0.0177694, 
  af[1, 3] -> 0.0485007, ag[1, 3] -> -0.0577328, af[2, 3] -> 0.069615,
   ag[2, 3] -> -0.113453, af[3, 3] -> 0.0786742, 
  ag[3, 3] -> -0.206362, af[4, 3] -> 0.0772869, ag[4, 3] -> -0.388605,
   af[5, 3] -> 0.0690158, ag[5, 3] -> -0.818868, 
  af[6, 3] -> 0.0575534, ag[6, 3] -> -2.14983, af[7, 3] -> 0.0458242, 
  ag[7, 3] -> -9.00525}}  

We can compare s1=.666052 with exact solution maxBid = 2/3. It looks like we have error of $6\times 10^{-4}$. Finally we plot numerical solution (red points) with exact solution (solid lines)

fe[s_] := 2 s/(1 + 3/4 s^2); ge[s_] := 2 s/(1 - 3/4 s^2);
lst1 = Table[{s, f0[s] /. solM1[[2]]}, {s, 0, .666, .01}]; lst2 = 
 Table[{s, g01[s] /. solM1[[2]]}, {s, 0, .666, .01}];

{Show[Plot[fe[s], {s, 0, maxBid}, AxesLabel -> {"s", "f"}], 
  ListPlot[lst1, PlotStyle -> Red]], 
 Show[Plot[ge[s], {s, 0, maxBid}, AxesLabel -> {"s", "g"}], 
  ListPlot[lst2, PlotStyle -> Red]]}