Re: how to plot nminimized result
- To: mathgroup at smc.vnet.net
- Subject: [mg113724] Re: how to plot nminimized result
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Tue, 9 Nov 2010 03:52:55 -0500 (EST)
tarun dutta wrote:
> On Nov 1, 2:59 pm, Daniel Lichtblau <d... at wolfram.com> wrote:
>> ----- Original Message -----
>>> From: "tarun dutta" <tarundut... at gmail.com>
>>> To: mathgr... at smc.vnet.net
>>> Sent: Sunday, October 31, 2010 2:09:22 AM
>>> Subject: how toplotnminimizedresult
>>> p = 0.01;
>>> q = 1;
>>> n = 5;
>>> CompNum[i_] := a[i] + I b[i];
>>> TCompNum[i_] := a[i] - I b[i];
>>> f = Sum[SetPrecision[
>>> Sqrt[i + 1] p CompNum[i] TCompNum[i + 1] +
>>> p CompNum[i + 1] TCompNum[i] + (q*i + i (i - 1)) CompNum[
>>> i] TCompNum[i], Infinity], {i, 0, n}];
>>> c = \!\(
>>> \*UnderoverscriptBox[\(\[Sum]\), \(i = 0\), \(n\)]
>>> \*SuperscriptBox[\(Abs\ [CompNum[i]]\), \(2\)]\) ==
>>> 1; c // TraditionalForm;
>>> dp = {a[n + 1] -> 0, b[n + 1] -> 0};
>>> var = Table[CompNum[i], {i, 0, n}] /. a_ + I b_ -> {a, b} // Flatten;
>>> prob = Join[ComplexExpand[f] /. a_ + I b_ -> {a^2 + b^2} /. dp, {c}];
>>> {val, res} =
>>> NMinimize[prob, var, MaxIterations -> 10000,
>>> AccuracyGoal -> 30]; // AbsoluteTiming
>>> this is my main program.for fixed value of p and q.I got the minimized
>>> value and also the variable like a[o],a[1]..etc for example(from above
>>> program)
>>> In[12]:= val
>>> Out[12]= 5.25767*10^-6
>>> In[13]:= res
>>> Out[13]= {a[0] -> -0.805443, b[0] -> 0.591251, a[1] -> -0.00824279,
>>> b[1] -> 0.0402093, a[2] -> 0.0000547477, b[2] -> -0.000202243,
>>> a[3] -> -1.11872*10^-7, b[3] -> 3.84618*10^-7,
>>> a[4] -> 6.22738*10^-10, b[4] -> 8.88902*10^-11,
>>> a[5] -> 3.02522*10^-10, b[5] -> 3.02899*10^-10}
>>> now I want to vary the value of p from 0 to 2 for a fixed value of q.q
>>> will also vary from 0 to 5.so, program will start taking the first
>>> value of q as 0 and scan p from 0 to 2 in steps 0.01.every time i will
>>> get the corresponding {val,res}.for example
>>> q=1 p=0.01 val 5.25767*10^-6 res a[0] -> -0.805443,
>>> b[0] -> 0.591251, a[1] -> -0.00824279,
>>> b[1] -> 0.0402093, a[2] -> 0.0000547477, b[2] ->
>>> -0.000202243,
>>> a[3] -> -1.11872*10^-7, b[3] -> 3.84618*10^-7,
>>> a[4] -> 6.22738*10^-10, b[4] -> 8.88902*10^-11,
>>> a[5] -> 3.02522*10^-10, b[5] -> 3.02899*10^-10}
>>> q=1 p=0.02 val= res =
>>> now i want to check theresultof res----if any a[i] and b[i] of all
>>> a[i] and b[i] is nearly equal to 1(or > o.1) and other a[i] are zero
>>> then print 'true' and also print the corresponding value of p and q
>>> if more than one a[i] and b[i] have value nearly equal to 1 then print
>>> 'false'
>>> from above example..
>>> q=1 p=0.01 res=True because only a[0] and b[0] have nearly = to=
> 1.we
>>> ignore all other a[i] and b[i] cause they have value in order of 10^-2
>>> or more.
>>> so there will be some kind of table as
>>> q=1 p=0.01 res true
>>> q=1 p=0.02 res true
>>> q=1 p=0.03 res true
>>> q=1 p=0.04 res false
>>> now we only consider only the last value of p for which we get
>>> res=true after that point we get false
>>> from above example we note the value of p=0.03 for q=1
>>> similarly
>>> q=1.1 res true p=0.01
>>> q=1.1 res true p=0.02
>>> q=1.1 res false p=0.03
>>> here we note the value p=0.02,q=1.1
>>> in this way we get a table of true value like
>>> q=1 p=0.03
>>> q=1.1 p=0.02
>>> ....
>>> .....
>>> now I want toplot(contour) q vs p.......
>>> this is my problem....so how will I do all this in mathematica...
>>> help
>>> if the problem is not clear to you people just mail me...
>>> with regards,
>>> tarun
>> If it helps anyone, below is code I sent to the original poster in privat=
> e email several days ago. It makes a function of the parameters (p,q) and a=
> lso uses a faster minimization (FindMinimum with interior point).
>> n = 5;
>> x[6] = 0;
>> x[i_] = re[i] + I*im[i];
>> conj[a_] := ComplexExpand[Conjugate[a]]
>> f[p_, q_] =
>> Together[Sum[
>> Sqrt[i + 1]*(p x[i]*conj[x[i + 1]] +
>> p*x[i + 1] conj[x[i]]) + (q*i + i (i - 1)) x[i]*
>> conj[x[i]], {i, 0, n}]];
>> c = Expand[Sum[x[i]*conj[x[i]], {i, 0, n}]] == 1;
>> v = Join[Array[re, n + 1, 0], Array[im, n + 1, 0]];
>> fmin[p_?NumericQ, q_?NumericQ] :=
>> FindMinimum[{f[p, q], c}, v, Method -> "InteriorPoint"]
>>
>> Even using FindMinimum instead of NMinimize, this can take time. With inc=
> rements of 1/10 instead of 1/1000, we have
>> In[72]:= Timing[
>> pts = Flatten[
>> Table[{p, q, fmin[p, q]}, {p, 0, 2, 1/10}, {q, 0, 5, 1/10}], 1];]
>> Out[72]= {75.598, Null}
>>
>> So a run over the full requested grid might be over two hours. Possibly u=
> sing sensible explicit initial values for the variables would help; I have =
> not tried that.
>> One probably would do better to use explicit loops on {q,p}, and terminat=
> e the inner loop whenever the condition above goes to False. I had not been=
> addressing that issue, hence the full table.
>> Daniel Lichtblau
>> Wolfram Research
>
> thanks daniel for replying once again....
> I need a global minima that is why I am using 'nminimize' instead of
> 'find minima'.the solution according to my problem which you cited
> here is not working properly.kindly once again you go through my
> problem...in your program I can not check the result..
> and with that I can not plot it also....
> plz help..
>
I doubt NMinimize is doing better than FindMinimum, for this range of
values. If and when I see examples where it is, I can reconsider that
aspect.
As for finding p-q pairs, this is straightforward programming.
findp[q_,eps_] := Module[
{p=0, cond=True, vals},
While[cond,
p += .01;
vals = v /. fmin[p,q][[2]];
If [Count[vals, a_ /; a>=eps] >=3, cond = False];
];
p - .01
]
One might speed this with a binary search but I'll skip that aspect.
The rest of the code I used is as I indicated above.
With a threshold (eps parameter) of .01, I get results as below. Here
the threshold is simply the value above which we consider variables to
be sufficiently "nonzero". When we have three such, we know p is too
large. This at least is my interpretation of what was indicated is needed.
In[11]:= qptab = Table[{q, findp[q, .01]}, {q, 1, 5, .1}]
Out[11]= {{1., 0.2}, {1.1, 0.21}, {1.2, 0.23}, {1.3, 0.24}, {1.4,
0.26}, {1.5, 0.27}, {1.6, 0.29}, {1.7, 0.3}, {1.8, 0.32}, {1.9,
0.33}, {2., 0.35}, {2.1, 0.36}, {2.2, 0.38}, {2.3, 0.39}, {2.4,
0.4}, {2.5, 0.42}, {2.6, 0.43}, {2.7, 0.45}, {2.8, 0.46}, {2.9,
0.48}, {3., 0.49}, {3.1, 0.51}, {3.2, 0.52}, {3.3, 0.53}, {3.4,
0.55}, {3.5, 0.56}, {3.6, 0.58}, {3.7, 0.59}, {3.8, 0.61}, {3.9,
0.62}, {4., 0.63}, {4.1, 0.65}, {4.2, 0.66}, {4.3, 0.68}, {4.4,
0.69}, {4.5, 0.71}, {4.6, 0.72}, {4.7, 0.73}, {4.8, 0.75}, {4.9,
0.76}, {5., 0.78}}
To get an idea of what the q vs p curve looks like, do
ListPlot[qptab]
Daniel Lichtblau
Wolfram Research