# 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[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?

• You might be interested in ParametricNDSolve[]. Mar 19 '18 at 22:55
• Does it help to solve these ODEs symbolically instead of numerically? Mar 20 '18 at 2:03
• @bbgodfrey I'd like to have a general solution for any ODE, even the ones not symbolically solvable. 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. Mar 21 '18 at 4:28
• By the way, what is the value of b? 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[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.

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