3
$\begingroup$

I have a numerical solution and I want to find the maximum value of the solution.

a = Rationalize[0.33390683175743086];
b = -0.15;
d = 5;
a1 = 
  u[y_] := 
    (Tanh[a*(y - d)] + Tanh[a*d])/(1 + Tanh[a*d]) + 
      b*Sqrt[3]*y/d*Exp[-1.5*(y/d)^2 + 0.5];
α = 0.04;
ω = Rationalize[0.0060230765816024628] + Rationalize[0.0097304407482613764]*I;
sol = 
  NDSolve[
    {(u[y] - ω/α)*(ϕ''[y] - α^2*ϕ[y]) - u''[y]*ϕ[y] == 0, 
     ϕ[80] == 1, ϕ'[80] == -I*α}, ϕ, {y, 0, 80}][[1, 1, 2]]
ϕ[y_] := sol[y]

So ϕ is the solution of my problem. I difined a new function called uu which is the negative of the derivative of ϕ:

uu[y_] := -ϕ'[y];

When I plot Abs[uu[y]] I have this:

enter image description here

Now the question is how to determine the maximum value of |uu|? Looking the figure it seems that the value is around for |uu| = 2.6.

Improve this question
$\endgroup$

4 Answers 4

Reset to default
7
$\begingroup$

You need to use arbitrary precision rather than machine precision in NDSolve by setting the WorkingPrecision. This then requires that all numeric values have at least that precision. Consequently, Rationalize or set to exact numbers all of your constants. To Rationalize long decimals you need to use a smaller tolerance. 0 will always work. 

a = Rationalize[0.33390683175743086, 0];
b = -3/20;
d = 5;
a1 = u[y_] := (Tanh[a*(y - d)] + Tanh[a*d])/(1 + Tanh[a*d]) + 
    b*Sqrt[3]*y/d*Exp[-3/2*(y/d)^2 + 1/2];
α = 1/25;
ω = Rationalize[0.0060230765816024628 + 0.0097304407482613764*I, 0];

Clear[ϕ];

sol = NDSolve[{(u[y] - ω/α)*(ϕ''[
           y] - α^2*ϕ[y]) - u''[y]*ϕ[y] == 0, ϕ[80] == 
      1, ϕ'[80] == -I*α}, ϕ, {y, 0, 80}, 
    WorkingPrecision -> 15][[1, 1, 2]];

ϕ[y_] := sol[y]

uu[y_] := -ϕ'[y]

Plot[Abs[uu[y]], {y, 0, 80}, PlotRange -> All]

enter image description here

Maximize[{Abs[uu[y]], 0 < y < 80}, y]

(* {2.63046649575507, {y -> 5.00764789850856}} *)

EDIT: Alternatively, rather than setting the WorkingPrecision in NDSolve; instead of using Maximize use NMaximize and set its WorkingPrecision.

NMaximize[{Abs[uu[y]], 0 < y < 80}, y, WorkingPrecision -> 15]

(* {2.63046694880644, {y -> 5.00764953131638}} *)

There is a slight difference in the result.

Also note that just setting the WorkingPrecision in Maximize does not work.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Thanks @BobHanlon. Solved. If you could explain me: Why I need to use arbitrary precision rather than machine precision in NDSolve? $\endgroup$
    – Mateus
    Commented Nov 10, 2017 at 20:41
  • 2
    $\begingroup$ Remove the option for WorkingPrecision and reevaluate; you will get an error message: "NMaximize::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations." Such problems are frequently due to loss of precision using machine precision. Arbitrary precision tracks and controls precision during calculations. The trade-off is that arbitrary precision is slower. Consequently, it is often only used if necessary. $\endgroup$
    – Bob Hanlon
    Commented Nov 10, 2017 at 20:53
5
$\begingroup$

Another answer copied from Daniel Lichtblau's. I updated it to use WhenEvent instead of EventLocator. The idea is to integrate the derivative of the function and use WhenEvent to detect extrema.

f[y_?NumericQ] = Abs[uu[y]]
vals = Reap[
    soln = y[x] /. First[
      NDSolve[{
        y'[x] == Evaluate[f'[x]],
        y[10] == (f[10]), 
        WhenEvent[y'[x] == 0, Sow[{x, y[x]}]]}, y[x], {x, .5, 30}]]][[2, 1]]

Then plot the function with the points stored in vals:

Plot[f[x], {x, 0, 30}, Epilog -> {PointSize[Medium], Red, Point[vals]}]

enter image description here

It then suffices to select the point of vals with the largest second value (Last[SortBy[vals, Last]]).

Improve this answer
$\endgroup$
3
$\begingroup$

One can extract the interpolated values from the InterpolatingFunction object, so here is another possibility. I've rewritten your ODE slightly:

sol = NDSolveValue[
    {(u[y] - ω/α)*(ϕ''[y] - α^2*ϕ[y]) - u''[y]*ϕ[y] == 0, ϕ[80] == 1, ϕ'[80] == -I*α},
    ϕ,
    {y, 0, 80}
];
sol //OutputForm

InterpolatingFunction[{{0., 80.}}, <>]

Then, the approximate maximum can be found from:

Max @ Abs @ (sol')["ValuesOnGrid"]

2.63037

This value will be slightly lower than the real maximum because the maximum probably occurs between grid points.

Improve this answer
$\endgroup$
2
  • $\begingroup$ That could be well combined with Henrik's answer :) $\endgroup$
    – anderstood
    Commented Nov 10, 2017 at 23:23
  • $\begingroup$ One could also consider using PeakDetect[] on the results of "Grid" and "ValuesOnGrid". $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Nov 11, 2017 at 0:58
2
$\begingroup$ 
vv[y_] := ϕ'[y] Conjugate[ϕ'[y]];
min = FindMaximum[{vv[t], 0. <= t <= 80}, {t, 2.6}]
Sqrt[min[[1]]]

{6.91936, {t -> 5.00765}}

2.63047

Improve this answer
$\endgroup$
1
  • $\begingroup$ But in this case I have to give an aproximation for "t" and the function FindMaximum gives then a local maximum. Is there a way to catch the global maximum without give an aproximation for 't'? $\endgroup$
    – Mateus
    Commented Nov 10, 2017 at 20:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.