Trying to generate a function from InterpolatingFunction with FunctionInterpolation

Here is the code (I’m sorry to give you such a complex code sample but I failed to reproduce the error with a simpler InterpolatingFunction):

(*constants*)
R0 = 8.314; M = 28.8 10^(-3); R = R0/M;
te0 = 298; pa0 = 101325; ρg0 = Pi 25 10^-4 pa0/(R te0);
ρs0 = 0.17; t0 = 0.005; x1 = 0.01; x2 = -10^-3;
Q = 3000000; A = 10^15; ea = 9 10^4; cg = (5/2) R; cs = 1.46 10^3;
kg = 2.61 10^-2; ks = 0.36; ac = 4/0.01; hc = 1500;

c1 = NDSolve[{((ArcTan[-10^5 x]/Pi) + 1/2) (ρs0 cs D[te[t, x], t] –
(0.005^2 Pi) ks D[te[t, x], x, x] – ρs0 Q A Exp[-ea/(R0 te[t, x])])
+ ((ArcTan[1000000 x]/Pi) + 1/2) (ρg0) cg D[te[t, x], t]
(+ρg0) (A ρs0/ρg0 Exp[-ea/(R0 te0)] x2 –
A ρs0/ρg0 Exp[-ea/(R0 te0)] x2/x1 x) (1 – Exp[-10^4 t])
cg D[te[t, x], x] (+ρg0) R te[t, x] (1 – Exp[-10^4 t])
(-A ρs0/ρg0 Exp[-ea/(R0 te0)] x2/x1) –
((0.005^2 Pi) kg D[te[t, x], x, x]) + (0.005^2 Pi) ac hc (te[t, x] – te0) == 0,
te[0, x] == te0, te[t, x2] == te0, te[t, x1] == te0},
{te[t, x]}, {t, 0, t0}, {x, x2, x1}, PrecisionGoal -> 3];

f[t_, x_] = te[t, x] /. c1; teb1[t_] = f[t, 0];
Plot[teb1[t], {t, 0, t0}]
Subscript[teb, i] =
FunctionInterpolation[
If[(D[(te[t, x] /. c1), x] /. x -> 0) < 100, ((te[t, x] /. c1) /. x -> 0) + 0.01,
((te[t, x] /. c1) /. x -> 0)],
{t, 0, t0}]

After I ran the code, the warning message generated…”The function did not evaluate to a real number”, why? Who can tell me what is wrong?

=================

  

 

Looks like condition …<100 in If is never satisfied. Is this correct? – Vitaliy Kaurov Aug 5 '12 at 7:36      Well…I thought I've seen an answer and I was already trying to comment it but it disappeared after I refresh this page…in fact <100 is already a simplified condition, the original one is more complex, but I think it doesn't matter, right? for D[(te[t, x] /. c1), x] /. x -> 0) is always a real number, if it’s not <100 it will >=100 and the value will turn to ((te[t, x] /. c1) /. x -> 0).
– xzczd
Aug 5 ’12 at 7:44

  

 

“is always a real number” – it’s not, but a real number in List; – that was the bug as @belisarius answered
– Vitaliy Kaurov
Aug 5 ’12 at 8:01

  

 

Please try to tidy up the code yourself next time. You will get more answers and less complains.
– Dr. belisarius
Aug 5 ’12 at 8:12

  

 

Hehe, I’ve really tried, but…you can see the code I added after your answer.
– xzczd
Aug 5 ’12 at 8:17

=================

1 Answer
1

=================

Change your last statement by

tt = FunctionInterpolation[
If[First@(D[(te[t, x] /. c1), x] /. x -> 0) < 100, ((te[t, x] /. c1) /. x -> 0) + 0.01,
((te[t, x] /. c1) /. x -> 0)],
{t, 0, t0}]

Then

Plot[tt[t], {t, 0, t0}]

Edit

The problem is that since

c1 == {{ te[t,x] -> Interp ….}}

the following

D[(te[t, x] /. c1), x] /. x -> 0 /. t -> t0/2

gives you

{-1586.49}

which has Head List, and can’t be compared with 100.

Edit 2

All xSolve commands return lists. You should select only one of the elements returned to use it.

Compare for example:

c1 = Solve[x == 3, x];
(x /. c1) == 3 (* Incorrect *)
(x /. c1[[1]]) == 3 (* Correct *)

(*
{3} == 3
True
*)

1

 

+1 – super quick debugging 😉
– Vitaliy Kaurov
Aug 5 ’12 at 7:55

  

 

Oh…but, in fact when I was trying to simplify my question, I formed this code:shiyan = FunctionInterpolation[Sin[t x], {t, 0, 10}, {x, 0, 2}];test1 = FunctionInterpolation[ If[(D[(Sin[t x]), x] /. x -> 1) < (x D[(Sin[t x]), x] /. x -> 1), Sin[t] + 0.1, Sin[t]], {t, 0, 10}];test2 = FunctionInterpolation[ If[(D[(shiyan[t, x]), x] /. x -> 1) < (x D[(shiyan[t, x]), x] /. x -> 1), shiyan[t, 1] + 0.1, shiyan[t, 1]], {t, 0, 10}];Plot[test1[t], {t, 0, 10}] Plot[test2[t], {t, 0, 10}] This code doesn’t run into the error (though it runs into another small error), why?
– xzczd
Aug 5 ’12 at 8:11

  

 

@xzczd Because the problem comes from the NDSolve output (c1)
– Dr. belisarius
Aug 5 ’12 at 8:17

  

 

@belisarius Well, you mean there’s a difference between the InterpolatingFunctions from NDSolve and the one from FunctionInterpolation?
– xzczd
Aug 5 ’12 at 8:23

  

 

@xzczd The solve commands return lists. Try c1 = Solve[x == 3, x]. Do you see the {{ }} wrapped around the result? Now compare x /. c1 with x /. c1[[1]]. That is what is happening with your code.
– Dr. belisarius
Aug 5 ’12 at 8:29