1
$\begingroup$

I have a problem when I'm trying to use NMinimize in my RecurrenceTable, like in this example :

int[a_, gamma_] = Simplify[Integrate[(-a^5 + a)/(gamma - a^3), a]];

moyenne = {1.0702438969345758`, 1.0990804105051308`};

alphaa[A1_, A0_, gammatheta_] = A1/(A0*gammatheta);

alphab[A1_, A0_, B0_, gammatheta_, 
   gamma_] = (gamma*(1 - A0^3/B0^3) + 
     alphaa[A1, A0, gammatheta]^3*(A0^3/B0^3))^(1/3);

gamma[gammar_, gammatheta_] = gammar/gammatheta;

B1[A1_, A0_, B0_, gammatheta_, 
   gammar_] = (alphab[A1, A0, B0, gammatheta, 
       gamma[gammar, gammatheta]]*B0)/(alphaa[A1, A0, gammatheta]*A0)*
   A1;


Relation[gammatheta_, gammar_, A1_, A0_, B0_] = 
  Simplify[int[
     alphab[A1, A0, B0, gammatheta, gamma[gammar, gammatheta]], 
     gamma[gammar, gammatheta]] - 
    int[alphaa[A1, A0, gammatheta], gamma[gammar, gammatheta]]];


recurrencetableA1B1 = 
 RecurrenceTable[{A1r[1] == 50, B1r[1] == 70, 
   A1r[n + 1] == A1 /. 
    NMinimize[{Re[
          Relation[moyenne[[1]], moyenne[[2]], A1, A1r[n], 
           B1r[n]]]^2 + 
        Im[Relation[moyenne[[1]], moyenne[[2]], A1, A1r[n], 
           B1r[n]]]^2, A1 > 50}, A1][[2]], 
   B1r[n + 1] == 
    B1[( A1 /. 
       NMinimize[{Re[
             Relation[moyenne[[1]], moyenne[[2]], A1, A1r[n], 
              B1r[n]]]^2 + 
           Im[Relation[moyenne[[1]], moyenne[[2]], A1, A1r[n], 
              B1r[n]]]^2, A1 > 50}, A1][[2]]), A1r[n], B1r[n], 
     moyenne[[1]], moyenne[[2]]] }, {A1r, B1r}, {n, 1, 30}]

The error reads :

NMinimize::nnum: The function value 0.0272898 Im[-(206882./A1r[n]^3)+2 Sqrt[3] ArcTan[Power[<<2>>] Plus[<<2>>]]-2 Sqrt[3] ArcTan[Power[<<2>>] Plus[<<2>>]]+<<9>>+Log[1.01788 +0.942684 Power[<<2>>]+0.873041 Power[<<2>>]]]^2+0.0272898 Re[-(206882./A1r[n]^3)+<<12>>+Log[1.01788 +0.942684 Power[<<2>>]+0.873041 Power[<<2>>]]]^2 is not a number at {A1} = {50.0915}.

I assume this problem appears because that NMinimize does not "read" the value of A1r[1] and B1r[1]. How can I solve this problem please ?

$\endgroup$
  • 2
    $\begingroup$ You can use NArgMin[...] instead of A1 /. NMinimize[...][[2]], this also issues error messages, but at least gives the answer. $\endgroup$ – Alx Dec 10 '19 at 3:35
1
$\begingroup$

Alx in a comment above correctly noted that using NArgMin[...] instead of A1 /. NMinimize[...][[2]] allowed the code in the question to produce the desired results, although with numerous error messages. Here, we offer two alternative approaches, one of which is much faster, and both of which produce no error messages. Begin be rewriting the code in the question as follows:

int[a_, gamma_] = Simplify[Integrate[(-a^5 + a)/(gamma - a^3), a]];
moyenne = {1.0702438969345758`, 1.0990804105051308`};
alphaa[A1_, A0_, gammatheta_] = A1/(A0*gammatheta);
alphab[A1_, A0_, B0_, gammatheta_, gamma_] = (gamma*(1 - A0^3/B0^3) + 
    alphaa[A1, A0, gammatheta]^3*(A0^3/B0^3))^(1/3);
gam[gammar_, gammatheta_] = gammar/gammatheta;
B1[A1_, A0_, B0_, gammatheta_, gammar_] = (alphab[A1, A0, B0, gammatheta, 
     gam[gammar, gammatheta]]*B0)/(alphaa[A1, A0, gammatheta]*A0)*A1;

Relation[gammatheta_, gammar_, A1_, A0_, B0_] = Simplify[int[
    alphab[A1, A0, B0, gammatheta, gam[gammar, gammatheta]], gam[gammar, gammatheta]] - 
    int[alphaa[A1, A0, gammatheta], gam[gammar, gammatheta]]];

recurrencetableA1B1 = RecurrenceTable[{A1r[1] == 50, B1r[1] == 70,
    A1r[n + 1] == NArgMin[{Abs@Relation[moyenne[[1]], moyenne[[2]], A1, A1r[n], 
        B1r[n]], A1 > 50}, A1], 
    B1r[n + 1] == B1[NArgMin[{Abs@Relation[moyenne[[1]], moyenne[[2]], A1, A1r[n],
        B1r[n]], A1 > 50}, A1], A1r[n], B1r[n], moyenne[[1]], moyenne[[2]]]},
    {A1r, B1r}, {n, 30}] 

Note that two changes have been made here, in addition to the change suggested by Alx: First, the function gamma has been renamed gam so that its name is different from that of the variable gamma used elsewhere in the code, which is safer. Second, Abs@Relation[...] is used instead of summing the squares of the real and imaginary parts of Relation[...], which is faster.

Running this code produces the results,

(* {{50., 70.}, {53.2283, 75.1996}, {56.6556, 80.7919}, {60.2934, 86.8073}, 
    {64.1534, 93.2783}, {68.2481, 100.24}, {72.5905, 107.731}, {77.194, 115.791}, 
    {82.073, 124.466}, {87.242, 133.802}, {92.7166, 143.852}, {98.5127, 154.67}, 
    {104.647, 166.317}, {111.136, 178.858}, {117.999, 192.362}, {125.253, 206.905}, 
    {132.918, 222.569}, {141.013, 239.442}, {149.559, 257.618}, {158.575, 277.2}, 
    {168.084, 298.3}, {178.106, 321.037}, {188.664, 345.541}, {199.78, 371.952}, 
    {211.476, 400.42}, {223.775, 431.108}, {236.7, 464.194}, {250.273, 499.868}, 
    {264.517, 538.335}, {279.453, 579.819}} *)

along with several NArgMin::nnum error messages in about 20 seconds on my computer.

The error messages can be eliminated by using the following instead of the last block of code above.

fm[a0_?NumericQ, b0_?NumericQ] := 
    NArgMin[{Abs@Relation[moyenne[[1]], moyenne[[2]], x, a0, b0], x > 50}, x];
recurrencetableA1B1fm = RecurrenceTable[{A1r[1] == 50, B1r[1] == 70,
    A1r[n + 1] == fm[A1r[n], B1r[n]], 
    B1r[n + 1] == B1[fm[A1r[n], B1r[n]], A1r[n], B1r[n], moyenne[[1]], moyenne[[2]]]},
    {A1r, B1r}, {n, 30}]

which produces the same answers in about the same amount of time.

FindRoot provides a much faster computation, however. Use

fr[a0_?NumericQ, b0_?NumericQ] := (x /. FindRoot[
    Relation[moyenne[[1]], moyenne[[2]], x, a0, b0], {x, 50, 500}]) // Abs;
recurrencetableA1B1fr = RecurrenceTable[{A1r[1] == 50, B1r[1] == 70,
    A1r[n + 1] == fr[A1r[n], B1r[n]], 
    B1r[n + 1] == B1[fr[A1r[n], B1r[n]], A1r[n], B1r[n], moyenne[[1]], moyenne[[2]]]},
    {A1r, B1r}, {n, 30}] 

which produces the same results in about 0.04 seconds.

Addendum

It is natural to suppose that modifying the arguments of Relation to read, Relation[gammatheta_, gammar_, A1_?NumericQ, A0_?NumericQ, B0_?NumericQ] = ... would eliminate the error messages. It does not. It does, however, increase the runtime by 60% for computations here based on NArgMin.

$\endgroup$
  • $\begingroup$ Thanks it's clear neat and smart ! $\endgroup$ – J.A Dec 19 '19 at 12:45

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.