3
$\begingroup$

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!

addendum

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}]

from user293787's very helpful answer is evaluable !

Thanks again

$\endgroup$
5
  • $\begingroup$ Why don't you use ParametricNDSolve? $\endgroup$
    – user21
    Jul 11, 2022 at 16:30
  • $\begingroup$ @user21 I wasn't sure how to use ParametricNDSolve with function . $\endgroup$ Jul 12, 2022 at 6:11
  • $\begingroup$ @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! $\endgroup$ Jul 12, 2022 at 6:44
  • $\begingroup$ No, I don't know what the issue might be. $\endgroup$
    – user21
    Jul 12, 2022 at 7:10
  • $\begingroup$ @user21 What a pity, I had the hope NeumannValue (changes v12.2 to v12.3) might be the issue $\endgroup$ Jul 12, 2022 at 7:19

3 Answers 3

6
$\begingroup$

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}]

Figure 1

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}]  

Figure 2

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

Figure 3

Finally we can improve code proposed by Ulrich Neumann as follows

Needs["NDSolve`FEM`"]

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}]
  

Fugure 4

$\endgroup$
5
  • $\begingroup$ 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) $\endgroup$ Jul 12, 2022 at 13:56
  • $\begingroup$ @UlrichNeumann You are welcome! It looks like code with NestList not working in version 13.1 as well. $\endgroup$ Jul 12, 2022 at 15:04
  • $\begingroup$ I played around with the options of NIntegrate until NestList worked. Don't know what changed from v12.2 to v13.1 $\endgroup$ Jul 13, 2022 at 10:33
  • $\begingroup$ Very clever to force integration to FEM mesh. Thanks. $\endgroup$ Jul 13, 2022 at 12:01
  • $\begingroup$ Until now my nestlist version is fastest on v12.2. Improvement is only ne essary for newer version $\endgroup$ Jul 13, 2022 at 13:29
4
$\begingroup$

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}]

enter image description here

$\endgroup$
3
  • $\begingroup$ (+1) but code needs to be improved as well. See update to my answer. $\endgroup$ Jul 13, 2022 at 11:30
  • $\begingroup$ +1 And Btw, please feel free to move the green checkmark to @AlexTrounev's much more comprehensive answer (or to your own). $\endgroup$
    – user293787
    Jul 13, 2022 at 11:33
  • 1
    $\begingroup$ @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! $\endgroup$ Jul 13, 2022 at 12:04
3
$\begingroup$

What about this?

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.

$\endgroup$
7
  • $\begingroup$ Nice idea , thanks. But unfortunately Table[U[i], {i, 1, 3}] doesn't evaluate. $\endgroup$ Jul 11, 2022 at 14:19
  • $\begingroup$ Really? It does on my machine. Takes about 2 seconds and returns a list of three InterpolatingFunctions. $\endgroup$
    – user293787
    Jul 11, 2022 at 14:22
  • $\begingroup$ What's your Mathematica version? (mine: v12.2) $\endgroup$ Jul 11, 2022 at 14:23
  • $\begingroup$ V12.3. But nothing here should depend on recent changes I think? $\endgroup$
    – user293787
    Jul 11, 2022 at 14:24
  • $\begingroup$ 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. $\endgroup$ Jul 11, 2022 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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