# Problem with the integration of this piecewise function

I encounter a problem when I try to integrate the piecewise function defined in the following code:

 ClearAll["Global*"];

(* Parameters of the problem*)

M = 1; (* Mass scale*)
A = M/4 ;(* Background Amplitude*)
a = A; (* Perturbation Amplitude*)
L = 0.1/M; (* Background Length scale*)
Lb = 30*L; (* Boundary distance*)
l = L/20; (* Perturbation Length scale*)
d = L/3; (* Spatial translation of the perturbation *)
tmax = 10*L; (* Maximal time*)
\[Kappa] = 3*M; (*Velociy field stiffness*)
\[Delta]\[Chi]0[x_] =
a/(1 + ((x - d)/l)^2) +
a/(1 + ((x + d)/l)^2); (* Initial perturbation profile*)

(* Construct an initial characteristic with the correct asymptotical \
behaviour *)

f1[x_] = -\[Kappa]*
x; (* Desired profile close to x=0 to create the caustic *)
f2[x_] = 1; (* Asymptotical profile *)
v[x_] = Piecewise[{
{f1[x], -L < x < L},
{f2[x], x < -Lb/2 ||  x > Lb/2},
{Interpolation[{{{-Lb/2}, f2[-Lb/2], f2'[-Lb/2], f2''[-Lb/2],
f2'''[-Lb/2], f2''''[-Lb/2]}, {{-L}, f1[-L], f1'[-L]}}][
x], -Lb/2 < x < -L},
{Interpolation[{{{Lb/2}, f2[Lb/2], f2'[Lb/2]}, {{L}, f1[L],
f1'[L]}}][x], L < x < Lb/2}}];

(* Define the initial condition for the field *)

s[x_] = M^2*Sqrt[(2 - 2*(v[x])^2)/(3 - (v[x])^2)];
q[x_] = Integrate[s[x], x,
Assumptions -> -Lb < x < Lb]; (* Integrate the spatial derivative*)


The piecewise function is called v[x] and it is defined explicitly on the intervals: [-L,L] ; ]-inf, Lb/2] and [Lb/2, +inf[ and defined as an interpolation on the remaining intervals. Then I define s[x] by using this piecewise function v[x]. When I try to define a function q[x] as the integral of s[x] , I have the following error:

Integrate::ilim: Invalid integration variable or limit(s) in -2.99988.


I have read the documentation on this error but I don't manage to relate this to my problem. Can anybody help me please?

Thanks a lot.

• Do you think you could reduce this into a minimal-working-example? – AccidentalFourierTransform May 12 '18 at 18:29
• I think you mean Integrate[s[x],{ x, -Lb , Lb}] rather than Integrate[s[x], x, Assumptions -> -Lb < x < Lb] – andre314 May 12 '18 at 18:49
• Integrate[s[x],{x,-.1,.1}] and Integrate[s[x],{x,1.5,3}] and Integrate[s[x],{x,-3,-1.5}] appear to work with the rest of your code unchanged, but anything which strays into your Interpolation zones fails. Does that give you any ideas what to start looking at? Defining q[x_] and then using x as a variable of integration worries me. – Bill May 12 '18 at 18:50
• @AccidentalFourierTransform : if you ignore the definition of the parameters, the example is really short. Please read from (* Construct an initial characteristic with the correct asymptotical \ behaviour *) – Targamas May 12 '18 at 20:42
• @andre Actually I want to define the function q[x] as an undefined integral but with [-Lb,Lb] as domains – Targamas May 12 '18 at 20:44

You can use NDSolve to compute a numerical antiderivative:

q1 = NDSolveValue[{qq'[x] == s[x], qq[-Lb] == 0}, qq, {x, -Lb, Lb}]
(* Some extrapolation warnings that are due to how NDSolve rewrites Piecewise[].
They are unimportant. *)

ListLinePlot@q0 Since the interpolations are between just two points, one can use InterpolatingPolynomial for slightly better results. (The OP's Interpolation leads to some small imaginary parts, due to a square-root of slightly negative numbers. For some reason, we don't get that with the InterpolatingPolynomial. The negative numbers are due to an overshoot around x == 1.5 (try 1 - (v[-1.5 + 10^-6])^2).

v1[x_] = Piecewise[{{f1[x], -L < x < L},
{f2[x], x < -Lb/2 || x > Lb/2},
{InterpolatingPolynomial[
{{{-Lb/2}, f2[-Lb/2], f2'[-Lb/2], f2''[-Lb/2], f2'''[-Lb/2],
f2''''[-Lb/2]}, {{-L}, f1[-L], f1'[-L]}}, x],
-Lb/2 < x < -L},
{InterpolatingPolynomial[{{{Lb/2}, f2[Lb/2], f2'[Lb/2]}, {{L}, f1[L], f1'[L]}}, x],
L < x < Lb/2}}];

Block[{v = v1},
q2 = NDSolveValue[{qq'[x] == s[x], qq[-Lb] == 0}, qq, {x, -Lb, Lb}]
]


It seems, MMA can't do Integration of squared interpolation functions.

Generate a new interpolationfunction of the squared function values.

Plot[(v[x])^2, {x, -5, 5}, PlotRange -> All] s2 = Interpolation[
Table[{x, M^2*Sqrt[(2 - 2*(v[x])^2)/(3 - (v[x])^2)]}, {x, -3, 3,
1/10}]];

Plot[s2[x], {x, -3, 3}, PlotRange -> All] is2[x_] = Integrate[s2[x], x]

Plot[is2[x], {x, -3, 3}]
` 