2
$\begingroup$

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?

$\endgroup$
1
  • 1
    $\begingroup$ 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[]. $\endgroup$ Feb 7 at 23:25
4
$\begingroup$

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

Figure 1

$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$ Feb 8 at 14:35
  • 2
    $\begingroup$ @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 $\endgroup$ Feb 8 at 16:32

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.