4
$\begingroup$

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?

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

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]

enter image description here

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]

enter image description here

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.

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

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.