# plotted function and Numerically integrate

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.

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!

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

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*


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}]


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}]


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}]


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.

• 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}] Dec 4 '17 at 14:41
• then Plot3D[f2sol[q, m5, [Theta]], {q, 0.8, 150}, {m5, 0, 150}] doesn't work. Dec 4 '17 at 14:41
• 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. Dec 5 '17 at 6:40