# Iterative pde solution NDSolve && NIntegrate

I try to solve a diffusion pde depending on a timedependent argument.

Step by step I have no problem to evaluate several solutions:

first step (timefunction 0):

tsim=10
U0 = NDSolveValue[{D[u[t, x], t] == D[u[t, x], {x, 2}] + NeumannValue[1, x == 0 ] - NeumannValue[u[t, x] , x == 1 ] + #  NeumannValue[1, x == 1 ], u[0, x] == 100 }, u , {t, 0, tsim}, {x, 0, 1} ] &  [0]


second step:

From result U0[t,x] a new timefunction int[t] is calculated

 int[t_?NumericQ] := NIntegrate[U0[t, x], {x, 0, 1}]


which is used to calculate the next solution.

third step (timefunction int[t]):

U1 = NDSolveValue[{D[u[t, x], t] == D[u[t, x], {x, 2}] + NeumannValue[1, x == 0 ] -
NeumannValue[u[t, x] , x == 1 ] + #  NeumannValue[1, x == 1 ], u[0, x] == 100 }, u , {t, 0, tsim}, {x, 0, 1} ] &  [int[t]]


Step 1-3 work fine ...

To easily automate these steps (iteration U[i][t,x]-> int[i][t]->U[i+1][t,x]->... I would like to use NestList but didn't succeed.

So far I tried

NestList[U =NDSolveValue[{D[u[t, x], t] ==D[u[t, x], {x, 2}] + NeumannValue[1, x == 0 ] -
NeumannValue[u[t, x] , x == 1 ] + #   NeumannValue[1, x == 1 ] , u[0, x] == 100 }, u , {t, 0, tsim}, {x, 0, 1} ] & ;
NIntegrate[U[t, x], {x, 0, 1}]
, 0, 3]


but Mathematica gives error messages NIntegrate::inumri: The integrand InterpolatingFunction[{{0.,10.},{0.,1.}},{5,4225,1,{106,41},{4,3},0,0,0,0,Indeterminate&,<<3>>},{{0.,<<9>>,<<96>>},<<1>>},{{{100.,100.,100.,100.,100.,100.,100.,100.,100.,100.,<<31>>},{169.706,29.1169,4.99567,0.857121,0.147059,0.0252312,0.00432868,0.000740845,0.000116391,-0.0000424981,<<31>>}},<<9>>,<<96>>},{Automatic,Automatic}] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0,1}}.

Thanks!

I would like to know from other Mathematica v12.2 users, wether the code

tsim=2; (* I assume this is a global variable with numerical value *)

U[i_/;i>=1]:=U[i]=With[{aux=int[i-1][t]},NDSolveValue[{D[u[t,x],t]==D[u[t,x],{x,2}]+NeumannValue[1,x==0]-NeumannValue[u[t,x],x==1]+aux* NeumannValue[1,x==1],u[0,x]==100},u,{t,0,tsim},{x,0,1}]];
int[i_][t_?NumericQ]:=If[i==0,0,NIntegrate[U[i][t,x],{x,0,1}]];

Table[U[i], {i, 1, 3}]


Thanks again

• Why don't you use ParametricNDSolve? Jul 11, 2022 at 16:30
• @user21 I wasn't sure how to use ParametricNDSolve with function . Jul 12, 2022 at 6:11
• @user21 Do you have an idea what could be the reason that user293787 's answer doesn't evaluate on Mathematica v12.2 (error NeumannValue...)? Thanks! Jul 12, 2022 at 6:44
• No, I don't know what the issue might be. Jul 12, 2022 at 7:10
• @user21 What a pity, I had the hope NeumannValue (changes v12.2 to v12.3) might be the issue Jul 12, 2022 at 7:19

I have v12.2 as well on my computer, therefore we can test code @ user293787, we have

tsim = 2;
U[i_ /; i >= 1] :=
U[i] = With[{aux = int[i - 1][t]},
NDSolveValue[{D[u[t, x], t] ==
D[u[t, x], {x, 2}] + NeumannValue[1, x == 0] -
NeumannValue[u[t, x], x == 1] + aux*NeumannValue[1, x == 1],
u[0, x] == 100}, u, {t, 0, tsim}, {x, 0, 1}]];
int[i_][t_?NumericQ] :=
If[i == 0, 0, NIntegrate[U[i][t, x], {x, 0, 1}]];


Visualization

Table[Plot3D[U[i][t, x], {t, 0, tsim}, {x, 0, 1},
PlotTheme -> "Marketing", ColorFunction -> Hue,
MeshStyle -> White], {i, 10}]


Code working nice and fast. We can compare solution with my code

tsim = 10;
int[0][t_] := 100; Do[
U1[i] = NDSolveValue[{D[u[t, x], t] ==
D[u[t, x], {x, 2}] + NeumannValue[1., x == 0.] -
NeumannValue[u[t, x], x == 1] +
int[i - 1][t] NeumannValue[1, x == 1], u[0, x] == 100},
u, {t, 0, tsim}, {x, 0, 1}];
int[i] = Interpolation[
Table[{t, NIntegrate[U1[i][t, x], {x, 0, 1}]}, {t, 0,
tsim, .1}]], {i, 5}]

Table[Plot3D[U1[i][t, x], {t, 0, tsim}, {x, 0, 1},
PlotTheme -> "Marketing", ColorFunction -> Hue,
MeshStyle -> White], {i, 5}]


Now we combine two solutions in one plot

tab = Table[{t, U[10][t, .5]}, {t, 0, 2, .1}]
Show[Plot[U1[5][t, .5], {t, 0, 2}], ListPlot[tab, PlotStyle -> Red]]


Finally we can improve code proposed by Ulrich Neumann as follows

Needs["NDSolveFEM"]

reg = ImplicitRegion[0 <= x <= 1, {x}]; mesH =
ToElementMesh[reg, MaxCellMeasure -> 2.5 10^-2];

tsim = 10; sol =
NDSolveValue[{D[u[t, x], t] ==
D[u[t, x], {x, 2}] + NeumannValue[1., x == 0.] -
NeumannValue[u[t, x], x == 1] +
100 NeumannValue[1, x == 1], u[0, x] == 100},
u[t, x], {t, 0, tsim}, Element[{x}, mesH]];

list = NestList[
NDSolveValue[{D[u[t, x], t] ==
D[u[t, x], {x, 2}] + NeumannValue[1., x == 0.] -
NeumannValue[u[t, x], x == 1] +
NIntegrate[# /. x -> x1, Element[{x1}, mesH]] NeumannValue[1,
x == 1], u[0, x] == 100}, u[t, x], {t, 0, tsim},
Element[{x}, mesH]] &, sol, 10]


Visualization

Table[Plot3D[list[[i]], {t, 0, tsim}, {x, 0, 1},
PlotTheme -> "Marketing", ColorFunction -> Hue,
MeshStyle -> White], {i, 10}]



• Thank you very much for your helpful answer and the confirmation "code runs on v12." . Now code also runs on my computer, no idea why... Interesting idea using Ìnterpoaltion( I tried FunctionInterpolation too) Jul 12, 2022 at 13:56
• @UlrichNeumann You are welcome! It looks like code with NestList not working in version 13.1 as well. Jul 12, 2022 at 15:04
• I played around with the options of NIntegrate until NestList worked. Don't know what changed from v12.2 to v13.1 Jul 13, 2022 at 10:33
• Very clever to force integration to FEM mesh. Thanks. Jul 13, 2022 at 12:01
• Until now my nestlist version is fastest on v12.2. Improvement is only ne essary for newer version Jul 13, 2022 at 13:29

Meanwhile I found a smart solution using NestList which is worth to be presented I think:

list=NestList[NDSolveValue[{D[u[t, x], t] ==D[u[t, x], {x, 2}] + NeumannValue[1., x == 0.] -  NeumannValue[u[t, x], x == 1] +
NIntegrate[#[t, x], {x, 0, 1} ,Method -> {Automatic, "SymbolicProcessing" -> 0}, AccuracyGoal -> 5 ] NeumannValue[1, x == 1],
u[0, x] == 100}, u, {t, 0, tsim}, {x, 0, 1}] &, 0 &, 10] // Quiet

Plot3D[ Evaluate[Through[Rest[liste][t, x]]], {t, 0, tsim}, {x, 0, 1}]


• (+1) but code needs to be improved as well. See update to my answer. Jul 13, 2022 at 11:30
• +1 And Btw, please feel free to move the green checkmark to @AlexTrounev's much more comprehensive answer (or to your own). Jul 13, 2022 at 11:33
• @user293787 It is not easy for me to follow your advice because your fast answer was the starting point of this interesting discussion. Thanks again! Jul 13, 2022 at 12:04

tsim=2; (* I assume this is a global variable with numerical value *)

U[i_/;i>=1]:=U[i]=With[{aux=int[i-1][t]},
NDSolveValue[{D[u[t,x],t]==D[u[t,x],{x,2}]+NeumannValue[1,x==0]-NeumannValue[u[t,x],x==1]+aux* NeumannValue[1,x==1],u[0,x]==100},u,{t,0,tsim},{x,0,1}]];
int[i_][t_?NumericQ]:=If[i==0,0,NIntegrate[U[i][t,x],{x,0,1}]];


Edit. The above code and Table[U[i],{i,1,3}] runs without warning in Version 12.3. But OP seems to get warnings related to NeumannValue[...,x==1]. Here is an alternative version where both NeumannValue[...,x==1] are combined into one, it also runs on my machine and seems to give identical results:

tsim=2; (* I assume this is a global variable with numerical value *)

U[i_/;i>=1]:=U[i]=With[{aux=int[i-1][t]},
NDSolveValue[{D[u[t,x],t]==D[u[t,x],{x,2}]+NeumannValue[1,x==0]+NeumannValue[aux-u[t,x],x==1],u[0,x]==100},u,{t,0,tsim},{x,0,1}]];
int[i_][t_?NumericQ]:=If[i==0,0,NIntegrate[U[i][t,x],{x,0,1}]];


Note that NeumannValue is a relatively new command, and was updated in Version 12.2.

• Nice idea , thanks. But unfortunately Table[U[i], {i, 1, 3}] doesn't evaluate. Jul 11, 2022 at 14:19
• Really? It does on my machine. Takes about 2 seconds and returns a list of three InterpolatingFunctions. Jul 11, 2022 at 14:22
• What's your Mathematica version? (mine: v12.2) Jul 11, 2022 at 14:23
• V12.3. But nothing here should depend on recent changes I think? Jul 11, 2022 at 14:24
• Trying to evaluate Table[U[i], {i, 1, 3}] I get errorNDSolveValue::fembceval: The boundary condition NeumannValue[1,x==1] could not be evaluated numerically. Jul 11, 2022 at 14:27