1
$\begingroup$

I want to speed up drawing of a function. The answer is here Numerically integrate a plotted function. But only one problem: it doesn't work for function of several variables. I tried to generalize it, but unfortunately it doesn't work. Here is my example.

I start with solving an equation:

eqn = j - Sqrt[q^2 + qp^2 - 2*q*qp*Cos[\[Theta]]] - 
         Sqrt[qp^2 + (1/2)*(16*m5^2 + ma^2 + mp^2 - 
                  Sqrt[(-(16*m5^2) - ma^2 - mp^2)^2 - 
                      4*(ma^2*mp^2 - 16*m5^2*qp^2)])] == 0; 
With[{gensol = Solve[eqn, qp]}, 
     Block[{m = 5.5, M = 300, Nc = 3, 
         c = ScientificForm[-44687.3983417778], 
         b = ScientificForm[161593.81818181818], 
         k1 = ScientificForm[16.485010961790245], 
         k2 = ScientificForm[-13.131344420001051], ma, mp, j}, 
       {j = Sqrt[q^2 + (1/2)*(16*m5^2 + ma^2 + mp^2 + 
                     Sqrt[(-(16*m5^2) - ma^2 - mp^2)^2 - 
                         4*(ma^2*mp^2 - 16*m5^2*q^2)])], 
          ma = Sqrt[-2*(M^2 - 2*(3*k1 + k2)*
                     (Sqrt[(c + M^2 + 2*m5^2)/(2*(k1 + k2))] + 
                          m*(b/(2*(c + M^2 + 2*m5^2))))^2 - c + 
                   2*m5^2)], mp = Sqrt[(2*b*m)/
                (Sqrt[(c + M^2 + 2*m5^2)/(2*(k1 + k2))] + 
                   m*(b/(2*(c + M^2 + 2*m5^2))))]}; 
        sols = gensol]]; 
qpC12 = Compile[{{q, _Complex}, {m5, _Complex}, 
         {\[Theta], _Complex}}, Evaluate[qp /. sols[[2]]], 
       RuntimeOptions -> "EvaluateSymbolically" -> False]; 
qp32 = Re[qpC12[q, m5, \[Theta]]]; 

Then I do:

f1[q_, m5_, \[Theta]_] := (q^2 qp32^2 Sin[\[Theta]]^2 )/
                          ((2 (2 \[Pi])^2)  Sqrt[-(2 q Cos[\[Theta]] qp32) + qp32^2 + q^2]) 

f3[q, m5, \[Theta]] = NDSolveValue[{D[f2[q, m5, \[Theta]], \[Theta]] == f1[q, m5, \[Theta]], 
f2[0, 0, 0] == 0}, f2, {\[Theta], 0, Pi/2}]

An errors apear:

NDSolveValue::noout: No functions were specified for output from NDSolveValue.

ConstantArray::ilsmn: Single or list of non-negative machine-sized integers expected at 
position 2 of ConstantArray[{TemporaryVariable$143960,TemporaryVariable$143961},-\[Infinity]].

Is there any way to do what I want? Thank you!

$\endgroup$
  • 2
    $\begingroup$ For a start, contemplate: NumericQ[1] (* True *) NumericQ[ScientificForm[1]] (* False *) $\endgroup$ – John Doty Dec 3 '17 at 18:59
  • 1
    $\begingroup$ NDSolveValue wants numeric input. If your only solve variable isTheta, you need to specify numeric values for q and m5. $\endgroup$ – Bill Watts Dec 3 '17 at 21:17
  • $\begingroup$ @BillWatts numeric values for q∈[0,150] and m5∈[0,150]. How may I give NDSolveValue numeric input in such case? $\endgroup$ – illuminates Dec 3 '17 at 22:57
  • $\begingroup$ * I want to plot Plot3D[f3[q, m5, \[Theta]], {q, 0.8, 150}, {m5, 0, 150}] $\endgroup$ – illuminates Dec 3 '17 at 23:01
2
$\begingroup$

I've corrected some things that will get you farther, but I still can't get you all the way.

First part looks ok.

eqn = j - Sqrt[q^2 + qp^2 - 2*q*qp*Cos[\[Theta]]] - 
    Sqrt[qp^2 + (1/2)*(16*m5^2 + ma^2 + mp^2 - 
         Sqrt[(-(16*m5^2) - ma^2 - mp^2)^2 - 
           4*(ma^2*mp^2 - 16*m5^2*qp^2)])] == 0;

The next part, get rid of the ScientificForms. They are for output only, and with them, your constants don't get assigned.

With[{gensol = Solve[eqn, qp]}, 
  Block[{m = 5.5, M = 300, Nc = 3, c = -44687.3983417778, 
    b = 161593.81818181818, k1 = 16.485010961790245, 
    k2 = -13.131344420001051, ma, mp, 
    j}, {j = 
     Sqrt[q^2 + (1/2)*(16*m5^2 + ma^2 + mp^2 + 
          Sqrt[(-(16*m5^2) - ma^2 - mp^2)^2 - 
            4*(ma^2*mp^2 - 16*m5^2*q^2)])], 
    ma = Sqrt[-2*(M^2 - 
         2*(3*k1 + 
            k2)*(Sqrt[(c + M^2 + 2*m5^2)/(2*(k1 + k2))] + 
             m*(b/(2*(c + M^2 + 2*m5^2))))^2 - c + 2*m5^2)], 
    mp = Sqrt[(2*b*m)/(Sqrt[(c + M^2 + 2*m5^2)/(2*(k1 + k2))] + 
         m*(b/(2*(c + M^2 + 2*m5^2))))]};
   sols = gensol]];

Then

qpC12 = Compile[{{q, _Complex}, {m5, _Complex}, {\[Theta], _Complex}},
   Evaluate[qp /. sols[[2]]], 
  RuntimeOptions -> "EvaluateSymbolically" -> False]

Its a good idea to test things at this point.

qpC12[75, 75, .5]
*65.7399 + 36.0074 I*

And add arguments to:

qp32[q_, m5_, \[Theta]_] = Re[qpC12[q, m5, \[Theta]]];

An interesting thing is your problem has almost not dependence on the variable m5,as shown by the plot.

Plot[Evaluate[qp32[40, m5, .1]], {m5, 0, 150}]

enter image description here

Using other values for q and theta exhibits similar behavior. Modify definition of f1 to include arguments of qp32.

f1[q_, m5_, \[Theta]_] := (q^2 qp32[q, 
      m5, \[Theta]]^2 Sin[\[Theta]]^2)/((2 (2 \[Pi])^2) Sqrt[-(2 q Cos[\[Theta]] qp32[q, m5, \[Theta]]) + qp32[q, m5, \[Theta]]^2 + q^2])

f1 now works

Plot[Evaluate[f1[q, 30, .1]], {q, 0, 150}]

enter image description here

The next part I have only been able to work by assigning q and m5 prior to the NDSolveValue

q = 19;
m5 = 20;
f2sol = NDSolveValue[{D[f2[\[Theta]], \[Theta]] == f1[q, m5, \[Theta]], f2[0] == 0}, f2, {\[Theta], 0, Pi/2}]

Plot[f2sol[\[Theta]], {\[Theta], 0, \[Pi]/2}]

enter image description here

And some values of q and m5 don't work at all. You know what you want better than I do and maybe this can help you to debug further.

$\endgroup$
  • $\begingroup$ Thank you, do you know how it can generalize to Plot3D? If i do f2sol[q, m5] = NDSolveValue[{D[f2[q, m5, [Theta]], [Theta]] == f1[q, m5, [Theta]], D[f2[q, m5, [Theta]], [Theta]] == 0}, f2, {[Theta], 0, Pi/2}] $\endgroup$ – illuminates Dec 4 '17 at 14:41
  • $\begingroup$ then Plot3D[f2sol[q, m5, [Theta]], {q, 0.8, 150}, {m5, 0, 150}] doesn't work. $\endgroup$ – illuminates Dec 4 '17 at 14:41
  • $\begingroup$ It looks like the instability is coming from f1. Execute Plot[f1[140, 20, \[Theta]], {\[Theta], 0, \[Pi]/2}, Exclusions -> None] to see what I mean. You can also try Plot3D[Evaluate[f1[q, 20, \[Theta]]], {q, 0, 150}, {\[Theta], 0, \[Pi]/2}]. The high discontinuity in f1 is creating problems with the whole solution. $\endgroup$ – Bill Watts Dec 5 '17 at 6:40

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.