Optimizing a parameter in an ODE

Question

I've solved a system of two ODEs using NDSolve which look like this

$\qquad y''-ky'=c, \quad y(t=0)=y_0, \quad y'(t=0)=\sin(a)$ $\qquad x''-kx'=0, \quad x(t=0)=x_0, \quad x'(t=0)=\cos(a)$

Where $c,q,k,a,y_0,x_0$ are constants.

I used NDSolve and obtained the y-x plot via ParametricPlot and added a Manipulate command to change parameters $k, a$ as the following:

k = 0.2;
x0 = 0;
y0 = 0;
c = 2;
T = 2;
Manipulate[
  ParametricPlot[
    Evaluate[
      {x[t], y[t]} /. 
        (NDSolve[
           {y''[t] - c - k*y'[t] == 0, y[0] == y0, y'[0] == Sin[a], 
            x''[t] - k*x'[t] == 0, x[0] == x0, x'[0] == Cos[a]}, 
          {x, y}, {t, 0, T}])], 
    {t, 0, T}, 
    PlotRange -> All], 
  {k, 0, 0.5}, 
  {a, 0, Pi/2}]

Now I want to find the optimum value for $a$ so that $x(t(y=b))$ for any given $k$ is maximum.

To find $t(y=b)$, I tried using NSolve as in the following, but somehow it's giving weird answers:

Manipulate[
  NSolve[
    Evaluate[
      {x[t], y[t]} /. 
        (NDSolve[
           {y''[t] + g - k*y'[t] == 0, y[0] == y0, y'[0] == v0*Sin[theta]}, 
           {y}, {t, 0, T}])] == b, t], 
  {k, 0, 0.5}, 
  {a, 0, Pi/2}]

What am I doing wrong? How does optimization works in Mathematica?

Answer 1

As I noted in a comment, this problem can be solved using NMaximize. Begin by solving the ODEs themselves.

k = 0.2; x0 = 0; y0 = 0; c = 2; T = 2; b = 3;
s = ParametricNDSolve[{y''[t] - c - k*y'[t] == 0, y[0] == y0, 
    y'[0] == Sin[a], x''[t] - k*x'[t] == 0, x[0] == x0, 
    x'[0] == Cos[a]}, {x, y}, {t, 0, T}, {a}];

where all constants except b are as in the question. Since b is not specified in the question, we choose it to be 3. If a flavor for the shapes of the solutions is desired, use Plot3D.

Plot3D[{(x /. s)[a][t], (y /. s)[a][t]}, {a, 0, Pi/2}, {t, 0, T}, 
  AxesLabel -> {a, t}, LabelStyle -> Directive[Bold, Black, Medium], ImageSize -> Large]

where x is orange and y is blue. Now plot the two solutions with ContourPlot, applying the constraint that y[a][t] == b.

Show[ContourPlot[(x /. s)[a][t], {a, 0, Pi/2}, {t, 0, T}, ContourLabels -> All], 
     ContourPlot[(y /. s)[a][t] == b, {a, 0, Pi/2}, {t, 0, T}, ContourStyle -> Black],
     FrameLabel -> {a, t}, LabelStyle -> Directive[Bold, Black, Medium],
     ImageSize -> Large]

We see in this case that the maximum value of x along the black curve y[a][t] == b is about 1.9 and lies near a == 0 and t == 1.64. More exact values are obtained from

NMaximize[{(x /. s)[a][t], (y /. s)[a][t] == b, 0 < a < Pi/2, 0 < t < T}, {a, t}]

(* {1.93756, {a -> 0., t -> 1.63756}} *)

This approach readily generalizes to more complicated ODEs, provided that the topology of the solutions is not too complicated.