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

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?

• Do you mean to compute maxBid as well using all constraints? Then it is a typical optimization problem, not for NDSolve[], but raise for NMinimize[]. Feb 7 at 23:25

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]]}


• Thank you! This indeed does work. However I have no idea what's going on. Is there a reference you would recommend to understand this technique? Feb 8 at 14:35
• @MichaelCurry There are a lot of papers about this method. I can recommend for beginners: Ü. Lepik (2009) Haar wavelet method for solving stiff differential equations, Mathematical Modelling and Analysis, 14:4, 467-481. Paper available on journals.vgtu.lt/index.php/MMA/article/view/6570 Feb 8 at 16:32