$RecursionLimit::reclim: Recursion problem

Question

Please, find bellow my original question. Here is all my code:

vm = 5;
r0 = 127.25;
di = r0*2;
de = 273;
v = 0.000000464;
ρf = 958;
ρt = 7854;
pr = 2.88;
kf = 0.656;
kt = 58;
cf = 4186;
ct = 470;
re = (vm*(2*r0/1000)/v)
e = 0.00015;
f = 1.325 (Log[0.27*(e/(2*r0/1000)) + 5.74*(1/re)^0.9])^(-2)
tinf = 298;
tme = 373;
αf = kf/(ρf*cf)
c1 = (-f*ρf*vm^2/(4*r0/1000))
g = 9.81;
β = 1/((tme + tinf)/2)
var = 0.000018;
prar = 0.703;
α = 0.000025;
kar = 0.028;

ym[r_] := (((r0/1000) - r)*vm*(f/8)^(0.5))/v;
εm[r_] := v*(0.4*ym[r]/6)*(1 + (r/(r0/1000)))*(1 + 2*(r/(r0/1000))^2);

momento = NDSolve[{(1/r)*D[(v + εm[r])*r*D[u[r], r], r] == (1/ρf)*c1, u'[0.0000001] == 0, u[(r0/1000)] == 0},u[r], {r, -0.0000001, (r0/1000)}, MaxSteps -> 10000];
mom[r_] = u[r] /. First[momento];

hi = (kf/(2*r0/1000))*((re - 1000)*pr*(f/8))/(1 + 12.7*(f/8)^0.5*(pr^(2/3) - 1))
ra = (g*β*((tme + tinf)/2 - tinf)*(de/1000)^3)/(var*α);
he = (kar/(de/1000))*(0.6 + (0.387*ra^(1/6))/(1 + (0.559/prar)^(9/16))^(8/27))^2

Here starts my problem. Mathematica isn't able to solve systems of PDEs with more than one independent variable. But Mathematica is able to solve each one individually. So, I have to solve iteratively both the two equation varying a boundary condition (by the way, it is a partial elliptic equation). But, at the moment this is not the problem.

l = 200;
tempo = 100;
Clear[tmed]; 
Clear[tmedpred]; 
Clear[ttf]; 
Clear[med];
Clear[tpred]; 
tmed[x_, t_] := tinf; (*Initial*)
med[x_, t_] := tinf; (*Initial*)
tmedpred[x_, t_] := tinf; (*Initial*)
tpred[x_, t_] := tinf; (*Initial*)
ttf[x_, r_, t_] := tinf; (*Initial*)
subrel = 0.4;
eps = 10^-4;

(*While[erro\[GreaterEqual]eps,*)

enf1 = NDSolve[{(ρf*cf*(D[tf[x, r, t], t] + mom[r]*D[tf[x, r, t], x])) == ((1/r)*
   D[(r*(αf + εm[r])*D[tf[x, r, t], r]), r]),tf[x, r, 0] == (tme - tinf)*Exp[-1000*x] + tinf,tf[0, r, t] == tme, (D[tf[x, r, t], r] /. r -> 0.0000001) == 
 0, -kf*(D[tf[x, r, t], r] /. r -> (r0/1000)) == (1 - Exp[-1000*x])*hi*(tmed[x, t] - ttf[x, r0/1000, t])},tf[x, r, t], {x, 0, l}, {r, 0.0000001, (r0/1000)}, {t, 0, tempo}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 13}}];
tff[x_, r_, t_] = tf[x, r, t] /. First[enf1];

med[x_, t_] := (NIntegrate[(mom[r]*tff[x, r, t]*r), {r, 0.0000001, r0/1000}])*2/((2/((r0/1000)^2)*NIntegrate[(mom[r]*r), {r, 0.0000001, r0/1000}])*(r0/   1000)^2));

tpred[x_, t_] := tmedpred[x, t] + subrel*(med[x, t] - tmedpred[x, t]);
tmed[x_?NumericQ, t_?NumericQ] := tpred[x, t];

ent1 = NDSolve[{ρt*ct*(D[tt[x, r, t], t]) == (1/r)*(D[(r*kt*D[tt[x, r, t], r]), r]), tt[x, r, 0] == tinf, -kt*(D[tt[x, r, t], r] /. r -> (r0/1000)) == (1 - Exp[-1000*x])*hi*(tmed[x, t] - ttf[x, r0/1000, t]), -kt*(D[tt[x, r, t], r] /. r -> (de/2000)) == (1 - Exp[-1000*t])*he*(ttf[x, de/2000, t] - tinf)}, tt[x, r, t], {x, 0, l}, {r, (r0/1000), (de/2000)}, {t, 0, tempo}];

tmedpred = tmed

Here is the problem. "tmedpred" shall return to the structure of the first NDSolve (enf1) and the code keep the loop up to an error value. Realize that the While structure is not in the code. I'm doing a manual evaluation.

The error: $RecursionLimit::reclim: Recursion depth of 1024 exceeded. >> The box message: "A very large output was generated. Here is a sample of it: (<<1>>) + subrel (med[100, 50] - tmedpred[100, 50])

The problem starts with the input tmedpred = tmed

(original question) I have two function:

med[x_?NumericQ, 
   t_?NumericQ] := ((NIntegrate[(mom[r]*tff[x, r, t]*r), {r, 
       0, 
       r0/1000}])*(2/((2/((r0/1000)^2)*
          NIntegrate[(mom[r]*r), {r, 0, r0/1000}])*(r0/
           1000)^2)));

and,

tmed[x_, t_] := tmedpred[x, t] + subrel*(med[x, t] - tmedpred[x, t]);

mom[r] and tff[x,r,t] are NDSolve output.

This functions are inside a While code and at the end of the loop,

tmedpred = tmed (tmedpred become tmed and the While code should keep iterating up to an error value)

tmedpred[187,100] (it is a simple evaluation of the function)

At this moment the warning

"$RecursionLimit::reclim: Recursion depth of 1024 exceeded."

And a box with the message

"A very large output was generated. Here is a sample of it:"

(<<1>>) + subrel (med[187, 100] - tmedpred[187, 100])) + subrel (med[187, 100] - tmedpred[187, 100])

Could some one help me with this problem? I was wondering if it is possible to create a function that is not connected to the original function but keep its original behavior(a kind of numerical function).

Regards!

Answer 1

Updating mutually dependent functions at each iteration

It seems like one important aspect which leads to the circular dependency problem here (pointed out earlier) is that you want to update your function definitions for tmed, tmedpred, ttf, med, and tpred at each iteration.

Your code uses global rules to create and define your functions (e.g. f[x_] := ...). Updating a definition that is created as a global rule is going to be hard or problematic when it is used in the definition of another global-rule function used to redefine it.

I think a good alternative approach in this case is to create your functions as Function objects. This will make it easier to update them with a simple Set (e.g. f = Function[{x, t}, ...]; oldf = f).

However, you'll also need to make sure that when you create and/or update your definitions, the underlying functions that the current definition depends on are "captured" in their plain Function-object form, rather than as their symbolic names; otherwise, the circular dependency problem will remain.

In order to do this "capturing", I've provided a helper function (captureValues) at the end of this post.

The following listing is a revised version of your second section of code. It has a couple of typos corrected, and it uses the approach described above to store and update your function definitions for each iteration:

l = 200;
tempo = 100;
Clear[tmed];
Clear[tmedpred];
Clear[ttf];
Clear[med];
Clear[tpred];

(* Use Function objects for iteratively updated function definitions *)
tmed = tinf &;(*Initial*)
med = tinf &;(*Initial*)
tmedpred = tinf &;(*Initial*)
tpred = tinf &;(*Initial*)
ttf = tinf &;(*Initial*)
(* Each time one of these functions is updated with another Function \
object, it should use a "snapshot" of any Function object or normal \
value it depends on. This snapshotting can be done with a helper \
function. *)

subrel = 0.4;
eps = 10^-4;

(* NumericQ-guarded rule-form version of tmed *)
Clear[tmedNum];
(* Always delegates to tmed; tmed itself will be updated each time *)
tmedNum[x_?NumericQ, t_?NumericQ] := tmed[x, t];

(*While[erro\[GreaterEqual]eps,*)

enf1 = NDSolveValue[{(\[Rho]f*
       cf*(D[tf[x, r, t], t] + mom[r]*D[tf[x, r, t], x])) == ((1/r)*
       D[(r*(\[Alpha]f + \[CurlyEpsilon]m[r])*D[tf[x, r, t], r]), r]),
     tf[x, r, 0] == (tme - tinf)*Exp[-1000*x] + tinf, 
    tf[0, r, t] == tme, (D[tf[x, r, t], r] /. r -> 0.0000001) == 
     0, -kf*(D[tf[x, r, t], r] /. r -> (r0/1000)) == (1 - 
        Exp[-1000*x])*hi*(tmed[x, t] - ttf[x, r0/1000, t])}, 
   tf, {x, 0, l}, {r, 0.0000001, (r0/1000)}, {t, 0, tempo}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 13}}];

ttf = enf1;

med =
  Function[{x, t},
    (NIntegrate[(mom[r]*ttf[x, r, t]*r), {r, 0.0000001, 
        r0/1000}])*2/((2/((r0/1000)^2)*
          NIntegrate[(mom[r]*r), {r, 0.0000001, r0/1000}])*(r0/1000)^2)
    ] // captureValues;

tpred =
  Function[{x, t},
    tmedpred[x, t] + subrel*(med[x, t] - tmedpred[x, t])
    ] // captureValues;

tmed = tpred;

ent1 = NDSolve[{\[Rho]t*
      ct*(D[tt[x, r, t], t]) == (1/r)*(D[(r*kt*D[tt[x, r, t], r]), 
        r]), tt[x, r, 0] == 
     tinf, -kt*(D[tt[x, r, t], r] /. r -> (r0/1000)) == (1 - 
        Exp[-1000*x])*
      hi*(tmedNum[x, t] - 
        ttf[x, r0/1000, t]), -kt*(D[tt[x, r, t], r] /. 
        r -> (de/2000)) == (1 - Exp[-1000*t])*
      he*(ttf[x, de/2000, t] - tinf)}, 
   tt[x, r, t], {x, 0, l}, {r, (r0/1000), (de/2000)}, {t, 0, tempo}];

tmedpred = tmed;

With the function evaluations now made independent of the function symbols that were used to define them, no circularity occurs:

In[231]:= tmedpred[187, 100]

Out[231]= 328.

However, fixing the circularity problem reveals an error arising from the input of the NDSolve expression for ent1:

I'm no expert with PDEs, so you may have to pick it up from here :-)

Capturing Function objects to prevent circular dependency

Here is an example of a helper function captureValues that can be used to create "independent" Function objects from ones that depend on other Function objects (and "normal" values i.e. values with non-empty OwnValues).

ClearAll[evaluateSymbols, captureValues]

evaluateSymbols = # //. s_Symbol :> With[{ss = s}, ss /; True] &;

captureValues[f : HoldPattern@Function[_ | PatternSequence[, __]]] :=
  evaluateSymbols@f

captureValues[f : HoldPattern@Function[params : {___Symbol}, __]] :=
  Block[params, evaluateSymbols@f]

Some example tests:

testf1 = Total@{##} &;

testf2 = testf1 @@ #1/#2 &;

testf3 = Function[{x, y}, x + y + testf2[{x, y}, y]];

x = 3; (* to show that same-named parameters are preserved *)

captureValues@testf3

(* Out: Function[{x, y}, x + y + ((Total[{##1}] &) @@ #1/#2 &)[{x, y}, y]] *)

Edit:

A general note about the "value capturing" helper function above: It only protects the outermost Function's parameters from replacement, and thus doesn't work in all cases when there are in-scope ownvalues for parameter symbols of nested functions. Furthermore, it doesn't take any other scoping construct into account. A full solution for this capturing idea is doable I think, but non-trivial (and probably a good suggestion to file with WRI Support).

Answer 2

The following set of expressions, which is taken from the relevant part of your code and simplified for discussion purposes, causes infinite recursion.

tpred[x_, t_] := .6 tmedpred[x, t] + .4 med[x, t]; 
tmed[x_?NumericQ, t_?NumericQ] := tpred[x, t]; 
tmedpred = tmed

tpred is defined in terms of tmedpred. tmed is defined in terms of tpred. Then tmedpred is made a synonym for tmed, making tpred a function of itself. Do you not see the circularity of this?