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;
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;


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:

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


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).


  • 1
    $\begingroup$ Please format your code and post a complete working example to reproduce your problem. $\endgroup$
    – Yves Klett
    Aug 7, 2014 at 6:38
  • 2
    $\begingroup$ Post your whole code which produce error $\endgroup$ Aug 7, 2014 at 6:43
  • 2
    $\begingroup$ Could some one help me? @molekyla777 ? $\endgroup$ Aug 7, 2014 at 21:12
  • $\begingroup$ @Bill, the publication "Advanced Numerical Differential Equation Solving in Mathematica" says that in the case of PDEs, systems are possible only for one independent variable. I tried, but the output is a message like "it is not an ODE system". A plot of "momento" show me that it is ok. There is a discontinuity very closed to 0 and that is why I'm integrating from 0.0000001 to (r0/1000). My only problem (I think) is when I do tmedpred = tmed and the $RecursionLimit error. $\endgroup$ Aug 8, 2014 at 2:29

2 Answers 2


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;

(* 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 *)
(* Always delegates to tmed; tmed itself will be updated each time *)
tmedNum[x_?NumericQ, t_?NumericQ] := tmed[x, t];


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, 
          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 - 
      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:

Errors from ent1 NDSolve evaluation

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[, __]]] :=

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 *)


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


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).

  • $\begingroup$ Thanks a lot @billisphere, this is really useful for me. As you can realize, I'm not able to do an elegant code in Mathematica. In other hand, Mathematica has been saving me in my master degree. By the way, I'm an expert in the "NDSolve ibcinc" warning. :) $\endgroup$ Aug 9, 2014 at 0:02
  • 1
    $\begingroup$ @LuisFernandoMoura, no problem, glad to have been able to help. You might want to take a look at the edit I just made to introduce an updated helper function. My first attempt at the helper function was bloated and wrong in an important case that it was trying to cover (when a Function has a parameter with the same name as a previously set symbol). It still worked for your example because that case didn't come up. Anyway, fixed now. And best of luck in your master's! $\endgroup$
    – William
    Aug 9, 2014 at 6:42
  • $\begingroup$ Hi @billisphere! I still have a problem. Your code works very well, but... do you have any ideal why the error NDSolveValue::ndnum: Encountered non-numerical value for a derivative at x == 0.`. ? This error is related to the the function tmedNum[x,t] because if I set tmedNum=0 or any constant, there is no problem. $\endgroup$ Aug 13, 2014 at 18:07
  • $\begingroup$ @LuisFernandoMoura, I'm not sure unfortunately, since I don't encounter that error here. In the code above, tmedNum is permanently defined to always "hand off" to tmed, so I can only guess a problem would be with whatever current function that tmed is defined as during the iteration in which you see that error. Try enabling the Debugger (under the Evaluation menu) and enable "Break at Messages" and "Show Stack", then run your code again. When the message occurs and the stack window fills, inspect the value of tmed to verify its value is what you were expecting. $\endgroup$
    – William
    Aug 15, 2014 at 14:38
  • $\begingroup$ p @billisphere, the problem is... even though it was defined tmedNum[x_?NumericQ, t_NumericQ] := tmed[x, t], if you evaluate tmedNum doing NumericQ[tmedNum], the output is "False". What's the best way to solve this? I don't know. Maybe creating an interpolation function from a table, but this way sounds to be an alternative not so smart. $\endgroup$ Aug 17, 2014 at 5:38

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?

  • $\begingroup$ Yes @m_goldber. Question: Is it possible to create a function (tmed) that keep its original behavior but is not connected to the other functions (tmedpred and med)? I kind of numerical function? Sorry if the question sounds to be silly. $\endgroup$ Aug 8, 2014 at 5:43
  • 1
    $\begingroup$ @LuisFernandoMoura. I would like to help you more, but I am not sure that I can. It has be 50 years since I last worked with PDEs. What little I recall about solving them is very much out of date. I can follow your code well enough to see the circular definitions I pointe out, but I don't fully understand it. $\endgroup$
    – m_goldberg
    Aug 8, 2014 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.