Optimizing a parameter in an ODE

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 == y0, y' == Sin[a],
x''[t] - k*x'[t] == 0, x == x0, x' == 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 == y0, y' == 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?

• You might be interested in ParametricNDSolve[]. – J. M. will be back soon Mar 19 '18 at 22:55
• Does it help to solve these ODEs symbolically instead of numerically? – bbgodfrey Mar 20 '18 at 2:03
• @bbgodfrey I'd like to have a general solution for any ODE, even the ones not symbolically solvable. – Alireza Mar 20 '18 at 9:38
• I believe that you are asking the following: With all constants but a specified, find the value of a that maximizes x[t] subject to the constrain that y[t] == b. Is that correct? If so, try,NMaximize. – bbgodfrey Mar 21 '18 at 4:28
• By the way, what is the value of b? – bbgodfrey Mar 21 '18 at 4:33

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 == y0,
y' == Sin[a], x''[t] - k*x'[t] == 0, x == x0,
x' == 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.

• Works like a charm! Thanks. – Alireza Mar 21 '18 at 7:53