1

I'm trying to solve a basic calculus of variations problem numerically, by minimizing an integral over an Interpolating Function.

Basically this:

data = {1, x2, x3, 1};
f = Interpolation[data];
obj = Integrate[f[x]^2, {x, 1, 4}];
NMinimize[obj, {x2, x3}]

But I get an error

NMinimize::nnum: The function value [...] is not a number at {x2,x3} = {0.918621,0.716689}.

The function seems to be ok:

Block[{x2 = 0, x3 = 0.5}, f[2.5]]
0.15625

What am I doing wrong?

Edit:

Thank you for the answers, but I still don't understand what I'm doing wrong, only how I can change my code so Mathematica won't throw an error.

In particular:

  1. f[1.5]^2 /. {x2 -> 1, x3 -> 2} and ((f /. {x2 -> 1, x3 -> 2})[1.5])^2 both evaluate without error, so why do I need to square my data instead of the InterpolationFunction?
  2. Why does converting from an InterpolatingFunction to a piecewise polynomial fix the problem?

Also note that the code is just a MWE. The objective function would in general be a complicated function of f and its derivatives. For example $$obj = f \sqrt{1+f'[x]^2}.$$ Now, I suppose I could resample that and generate a new InterpolationFunction and then extract a piecewise polynomial of that. But why do I have to tie myself in knots like this? I can evaluate my objective at every point, so NMinimize should be able to do so as well, shouldn't it?

I'd like to treat the InterpolationFunction as I would any continuous function. I thought that was the point.

As I mentioned I'm trying to numerically approximate a continuous solution, so the solution needs to be practical for having n variables instead of just 2. The approximation would of course improve as n becomes larger.

I suspect this has something to do with evaluation order, but none of I'm doing violates what I know about the Wolfram Language.

Edit #2:

Here is the actual (non-MWE) problem I'm trying to solve:

npts = 10;
vars = Unique[Table[x, npts]];
data = {1, Splice[vars], 1};
f = Interpolation[data, InterpolationOrder -> 3];
f2 = FunctionInterpolation[
   f[x] Sqrt[1 + f'[x]^2], {x, 1, Length[data]}];
obj = Integrate[f2[x], {x, 1, Length[data]}];
cons = Flatten[Table[{0 < x, x < 1}, {x, vars}]];
min = NMinimize[{obj, Splice[cons]}, vars];
ListPlot[data /. min[[2]]]

The solution should approximate a Catenary.

user2910571
  • 13
  • 3
  • 1
    Your problem lies with trying to square the InterpolationFunction. Try something like data1 = data^2; f1 = Interpolation[data1]; obj1 = Integrate[f1[x], {x, 1, 4}]; NMinimize[obj1, {x2, x3}]. – bbgodfrey Aug 02 '21 at 12:51
  • If the data in the Q is your actual use-case, convert the interpolation to a piece wise polynomial. You can search on site for how to do that. For instance, https://mathematica.stackexchange.com/a/212753/4999 might do the trick. – Michael E2 Aug 02 '21 at 13:51

1 Answers1

1

Since you want numerical optimization, use NIntegrate and tell NMinimze x2,x3 are numerical values.

But be aware that it depends on the choosen InterpolationOrder 

data = {1, x2, x3, 1};
f = Interpolation[data, InterpolationOrder -> 3];
obj[x2_?NumericQ, x3_?NumericQ] := NIntegrate[f[x]^2, {x, 1, 4}];
nmin = NMinimize[obj[x2, x3], {x2, x3}]
Plot[f[x]^2 /. nmin[[2]], {x, 1, 4}]

(*   {0.5, {x2 -> -0.111111, x3 -> -0.111111}}   *)


data = {1, x2, x3, 1};
f = Interpolation[data, InterpolationOrder -> 1];
obj[x2_?NumericQ, x3_?NumericQ] := NIntegrate[f[x]^2, {x, 1, 4}];
nmin = NMinimize[obj[x2, x3], {x2, x3}]
Plot[f[x]^2 /. nmin[[2]], {x, 1, 4}]


(*   {0.6, {x2 -> -0.2, x3 -> -0.2}}   *)

Edit

to do things analytically with Integrate, generate a second interpolation function f2 of f[x]^2 with FunctionInterplation . But you may get a slight loss of precision, play with Options[FunctionInterpolation] to get best results, e.g. higher AccuracyGoal or more InterpolationPoints...

f2 = FunctionInterpolation[f[x]^2, {x, 1, 4}]//Quiet

Now you can integrate with Integrate and even minimize with Minimize 

int2 = Integrate[f2[x], {x, 1, 4}] // Simplify

(*   (3 (349985 + 771384 x2^2 + 150015 x3 + 771384 x3^2 - 
     3 x2 (-50005 + 64261 x3)))/2000000   *)


min = Minimize[int2, {x2, x3}]


(*   {997/1999, {x2 -> -(2009/17991), x3 -> -(2009/17991)}}   *)


Plot[f2[x] /. min[[2]], {x, 1, 4}]
Akku14
  • 17,287
  • 14
  • 32
  • The main problem is, that Integrate can only work with the direct interpolated function f[x], but not with the squared f[x]2^2. But NIntegrate can. But you can go another way to get fully analytical solutions. ........ See my edit. – Akku14 Aug 02 '21 at 20:19
  • Your use of Quiet silences errors which look pertinent. I wonder if the difficulties I'm having getting the correct solution are related. – user2910571 Aug 03 '21 at 00:00
  • Since you offer variables x2 and x3, FunctionInterpolation must complain, but it works nevertheless ( with my MMA version 8.0). The function values at interpolation points depend on x2,x2 as you see with f2[[4]] // Simplify . What do you mean with "correct" solution. As i said, it depends on interpolation order. If you want to get the same result with f2 and Integrate as with f and NIntegrate, choose high interpolation points, e.g. f2 = FunctionInterpolation[f[x]^2, {x, 1, 4}, InterpolationPoints -> 100] – Akku14 Aug 03 '21 at 03:59