[Date Index]
[Thread Index]
[Author Index]
Re: NMinimize a function of NMaximize
*To*: mathgroup at smc.vnet.net
*Subject*: [mg79990] Re: [mg79947] NMinimize a function of NMaximize
*From*: Daniel Lichtblau <danl at wolfram.com>
*Date*: Fri, 10 Aug 2007 01:48:37 -0400 (EDT)
*References*: <200708090924.FAA21421@smc.vnet.net>
sapsi wrote:
> I have a function
> g[x_]:=1 (Pi - 2 x*Sqrt[1 - x^2] - 2 ArcSin[x])
> d[x_, o_, n_] := Abs[g[x] - Normal[Series[g[y], {y, o, n}]]] /. y -> x
>
> d[] is the absolute difference of g[] and its Taylor series
> approximation
>
> dm[p_] := NMaximize[{d[x, p, 3], 0 <= x <= 1}, x][[1]]
>
> dm[] is the maximum absolute difference between the Taylor series
> approximation(3rd order) around p and g[].
> I wish to find that value of 'p' that minimizes this maximum
> difference - i tried plotting and can see where the minima occurs but
> would like the exact value. So, it thought this would work
>
> NMinimize[{NMaximize[{ d[x, p, 3], 0 <= x <= 1}, x], 0 <= p <= 1}, p].
> Instead i get errors (briefly)
>
> NMaximize::nnum: The function value -Abs[-2.41057-(2 (<<19>>-p)^2 \
> p)/Sqrt[1-p^2]+2 p Sqrt[1-p^2]+(2 (<<1>>)^3)/(3 Sqrt[1-<<1>>] \
> (-1+p^2))-(4 (0.652468-p) (-1+p^2))/Sqrt[1-p^2]+2 ArcSin[p]] is not a
> \
> number at {x} = {0.652468}. >>
> NMaximize::nnum: "The function value \
> -Abs[-2.41057-(2\(<<19>>-p)^2\p)/Sqrt[1-p^2]+2\ p\ Sqrt[1-p^2]+(2\(<<\
> 1>>)^3)/(3\Sqrt[1-<<1>>]\(-1+p^2))-(4\(0.652468-p)\(-1+p^2))/Sqrt[1-p^
> \
> 2]+2\ ArcSin[p]] is not a number at {x} = {0.6524678079740285`}."
> NMaximize::nnum: The function value -Abs[-2.41057-(2 (<<19>>-p)^2 \
> p)/Sqrt[1-p^2]+2 p Sqrt[1-p^2]+(2 (<<1>>)^3)/(3 Sqrt[1-<<1>>] \
> (-1+p^2))-(4 (0.652468-p) (-1+p^2))/Sqrt[1-p^2]+2 ArcSin[p]] is not a
> \
> number at {x} = {0.652468}. >>
> NMinimize::nnum:
>
> Can anyone provide any pointers on how to find the minimum of dm[]?
> Thank you for your time
> Saptarshi
>
Getting the code to handle the two-level optimization is relatively
easy. You could do, for example:
g[x_] := (Pi - 2*x*Sqrt[1 - x^2] - 2*ArcSin[x])
d[x_, o_, n_] := Abs[g[x] - Normal[Series[g[y], {y,o,n}]]] /. y->x
maxp[x_,p_Real,n_] := First[NMaximize[{d[x,p,3], 0<=x<=1}, x,
Method->DifferentialEvolution]]
Timing[{min,pbest} = NMinimize[{maxp[x,p,3], 0<=p<=1}, p]]
Problem is it will be quite slow, and moreover not give a reliable
result. The reason for the latter deficiency is that the inner
optimization will often go to the "wrong" side of the x interval, and
not find the correct max.
A simple way around this is to realize (by computations, if you like,
which is what I did) that for given p the max is always at one of the
endpoints, either x=0 or x=1. So you can save both time and mistakes by
checking those directly.
I recast so that some computations do not get repeated. In particular
I'll start with the generic Taylor series.
d[x_, o_] = Abs[2*o*Sqrt[1 - o^2] -
(4*(-1 + o^2)*(-o + x))/Sqrt[1 - o^2] -
(2*o*(-o + x)^2)/Sqrt[1 - o^2] + (2*(-o + x)^3)/
(3*Sqrt[1 - o^2]*(-1 + o^2)) - 2*x*Sqrt[1 - x^2] +
2*ArcSin[o] - 2*ArcSin[x]]
maxp[p_Real] := Max[d[0,p],d[1,p]]
In[13]:= Timing[{min,pbest} = NMinimize[{maxp[p], 0<=p<=1}, p]]
Out[13]= {1.6321, {0.0635301, {p -> 0.593819}}}
If you check d[...,p] at this p value you will find that the values are
equal at endpoints x = {0,1}. So we find the min p just where the
endpoints balance.
By the way, this is distantly related to a bilevel optimization used for
Trefethen's 2002 SIAM challenge problem #5. There is a zip file at
http://library.wolfram.com/infocenter/Conferences/5353/
It contains a Mathematica notebook with code to handle that optimization.
Daniel Lichtblau
Wolfram Research
Prev by Date:
**Re: Re: show residuals of circle fit**
Next by Date:
**Re: Version 6 "Mathematica Book" - updated and expanded**
Previous by thread:
**NMinimize a function of NMaximize**
Next by thread:
**Re: NMinimize a function of NMaximize**
| |