Optimization of ODE with respect to the initial condition

Question

One has a (system) of ODEs with a one-parameter family of initial conditions. For example,

f[kk_] := N[aa = kk; s = NDSolve[{y''[t] == y[t]^2, y[1] == aa, y'[1] == 1},
                                  y, {t, 2}]; Evaluate[(y[2] /. s)][[1]]]

where k is the parameter, and one sees $f(k)$ as a function of $k$ works fine, as one can try to plot

Plot[f[x], {x, -3, 0}]

and obtain

However, when one tries to find the min, using FindMinimum or NMinimize, it does not work as it seems Mathematica would first perform symbolic calculation with

NDSolve[{(y''[t] == -2y[t], y[1] == x, y'[1] == 1}, y, {t, 2}]

which results an error as the initial condition is not a number, and returns y[2.] before continuing evaluating NMinimize[y[2.], x].

Does anyone know how to get around it, or some better ways to do the optimization?

Answer 1

Use ParametricNDSolve or it's cousin ParametricNDSolveValue:

s = ParametricNDSolveValue[
    {y''[t] == y[t]^2, y[1] == aa, y'[1] == 1},
    y[2],
    {t, 2},
    aa
];

NMinimize[s[aa], aa]

{0.0964433, {aa -> -1.65429}}

Answer 2

Harnessing symbolic functionality one can get a glimpse of beauty of the mathematical world and a quantum of deeper understanding, therefore often it is of high importance and preferable to exploit symbolic rather than numeric tools, even though in certain cases the latter would yield more. Thus it is a reasonable practice to put forward exact methods and use approximate ones only if symbolic functionality fails.

We'll demonstrate that the function to be minimized is quite involved $$y(t,a)=6\;\wp\bigg(t-1+\wp^{-1}\Big(\frac{a}{6};\;0,\frac{a^3}{54}-\frac{1}{36}\Big);\;0,\frac{a^3}{54}-\frac{1}{36}\bigg)$$ where $\wp$ is the Weierstrass elliptic function and it can be done with NMinimize as well as with FindMinimum, however getting the actual result involves some insight into the matter because of suboptimal handling of elliptic functions and integrals in Mathematica (a marvelous trademark figured out by J.M.). Encountering this issue wherever below one will be warned with this sign $**$.

Let's start with a few quite natural observations. Having equations of the form y''[t] == y[t]^2 it is a good idea to multiply it by y'[t] and integrate it once using initial (boundary) conditions. $0=y''(t)-y(t)^2\; $ thus $\;0=y'(t) y''(t)-y(t)^2 y'(t)\;$ i.e. $\;0=1/2 \;y'(t)^2 -1/3\; y(t)^3 +c$ and using initial conditions we find: $\;c=-1/2 + a^3/3 $. Rewriting $\;y(t)=6 w(t)$ one obtains: $$w'(t)^2-4w(t)^3+ \frac{a^3}{54} - \frac{1}{36}=0$$ This is a standard form of a differential equation for the Weierstrass elliptic function $w'(t)^2-4w(t)^3+ g_2 w(t)+g_3=0\;$ with $\;g_2=0$ and $g_3= \frac{a^3}{54} - \frac{1}{36}$. Mathematica can $**$ solve this more general form on the symbolic level without specifying initial conditions:

w[t]/. DSolve[ w'[t]^2 - 4 w[t]^3 + g2 w[t] + g3 == 0, w[t], t]

{ WeierstrassP[ t - C[1], {g2, g3}], WeierstrassP[ t + C[1], {g2, g3}]}

however one might be misled using it $**$ with $g_2=0$.

Now, using once more the initial conditions we find

C[1] == 1 - InverseWeierstrassP[ a/6, {0, a^3/54 - 1/36}] and finally:

y[t_, a_] := 6 WeierstrassP[ t - 1 + InverseWeierstrassP[ a/6, {0, a^3/54 - 1/36}],
                             {0, a^3/54 - 1/36}]

when plotting this function $**$ one uses Re[ y[2,a]] rather than y[2,a] alone, since the latter produces $**$ some nonvanishing imaginary perturbations.

Plot[ Re @ y[2, a], {a, -3, 0}, AxesOrigin -> {0, 0}, PlotStyle -> Thick]

$**$ Minimize[{y[2, a], -3 < a < 0}, a] doesn't work, however $**$ NMinimize[{ Re @ y[2, a], -3 < a < 0}, a] works fine and alternatively an excellent approximation $**$ with FindMinimum as you liked can be found:

FindMinimum[ Re @ y[2, a], {a, 0}] // Quiet

{0.0964433, {a -> -1.65429}}   

If we'd like studying full complexity of the function y, a useful tool might be e.g.

Manipulate[ Plot[ Re @ y[t, a], {a, -10, 10}, AxesOrigin -> {0, 0}, PlotStyle -> Thick], 
           {t, -30, 30}]

which wouldn't be readily available with a numerical approach.

Another nice application and a bit more detailed explanation of a few issues related to elliptic functions one can find in e.g. Integrate yields complex value, while after variable transformation the result is real. Bug?. Elliptic functions are solutions to certain nonlinear ordinary differential equations and being implemented in Mathematica 1.0 they might be slightly revised, although some useful functionality like WeierstrassHalfPeriodW1 is implemented also in verion 11.2.

Answer 3

A very simple fix for your problem is to prevent f from evaluating without numeric arguments. While I'm at it, I'll also rewrite f in a slightly nicer way:

Clear[f];
f[kk_?NumericQ] := Module[{sol, y},
  sol = First[
    N[NDSolve[{y''[t] == y[t]^2, y[1] == kk, y'[1] == 1}, y, {t, 2}]]
  ];
  y[2] /. sol
];

The _?NumericQ pattern on f prevents the function from evaluating symbolically. In this form, NMinimize and FindMinimum work immediately:

NMinimize[f[x], x]
FindMinimum[f[x], {x, 0}]

Finally, I'd like to point out that your use of Evaluate does nothing in your code. Evaluate is used to override a Hold attribute. For example:

Hold[Range[3], Evaluate[Range[3]]]