2

I have problem with my PDE. The problem is as bellow:

enter image description here I used "WhenEvent" operator to solve it. The code is: 

a = 255
b = 2.5
heat1 = NDSolveValue[{D[u[t, x], 
  t] - (0.0000001 D[u[t, x], x, x] + b*Exp[a *x]) == 
NeumannValue[0, True] + NeumannValue[-7 (u[t, x] - 25), x == .05],
u[0, x] == 25, 
WhenEvent[
u[t, x] > 250, {a = 0.01*a, b = .2 b, "RestartIntegration"}]}, 
u, {t, 0, 600}, {x, 0, .05}, 
Method -> {"MethodOfLines", "TemporalVariable" -> t, 
"SpatialDiscretization" -> {"FiniteElement", 
  "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}}]
pa = Plot[Evaluate[heat1[t, 0.05]], {t, 0, 600}, PlotRange -> All, 
AxesLabel -> {t, "T(.05,t)"}]

After running this code, I got the following errors.

" NDSolveValue::nbnum1: The function value        InterpolatingFunction[{{0.,0.05}},{5,4225,0,{101},{3},0,0,0,0,Automatic,{},{},False},{<<1>>},{25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,<<15>>,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,25.,<<51>>},{Automatic}][x]>250 is not True or False when the arguments are {0.,{<<1>>},{2.50058,3.22634,4.16337,5.37266,6.93321,8.94703,11.5458,14.8994,<<36>>,186525.,240704.,310619.,400841.,517266.,667489.,<<51>>}}. >>
General::stop: Further output of NDSolveValue::nbnum1 will be suppressed during this calculation. >> "

How can I fix it? Thank you

user21
  • 39,710
  • 8
  • 110
  • 167
Hamed Peyrovedin
  • 25
  • 6
  • At a given time, u[t, x] is an array of values. Which do you wish to compare with 250? Perhaps, the largest value? – bbgodfrey Dec 14 '16 at 20:53
  • @bbgodfrey that you for your comment. As you said, at a given time, u[t, x], is an array of values. I want to check all of them and see if this value is higher than 250 or not. if it is lower that value, the solution is OK. But if it is higher than 250, it should be solved based on the corrections that mentioned. Now how can I do this? – Hamed Peyrovedin Dec 14 '16 at 21:10
  • 1
    I have fixed the coding problems by defining a and b as discrete variables, and making related changes. However, your coded boundary condition at x == 0.05 is not the same as your boundary condition written in LaTex and is causing severe difficulties. – bbgodfrey Dec 14 '16 at 21:50
  • @bbgodfrey Can you please send the code here? Which part of the boundary condition is not true? – Hamed Peyrovedin Dec 14 '16 at 22:08
  • 1
    The LaTex boundary condition involves a temporal derivative, but the coded boundary condition involves a spatial derivative. Moreover, the spatial derivative for x just less than .05 is positive and large, but the coded boundary condition at x == 0.5 has a negative and large spatial derivative. This incompatibility causes rapid variation in u near x == 0.5. – bbgodfrey Dec 14 '16 at 22:13
  • In the latex formulation f(x,u) is a discontinuous function of x at any t, because u is a function of x. Is that what you intend? Your WhenEvent appears to be attempting switch a,b for all x based on some critera. Can you clarify when you intend. – george2079 Dec 14 '16 at 22:22
  • @bbgodfrey That is true, it was my mistake. I edited the LaTex boundary condition. Is it OK now? – Hamed Peyrovedin Dec 14 '16 at 22:23
  • @george2079 The question is based on the latex formulation. I think I did a mistake in my code. consider a system that the governing equation for that system is given in latex formulation, the heterogeneous term is not the same for all u(x,t), if u become more than 250, the heterogeneous term will be changed according to latex formulation. – Hamed Peyrovedin Dec 14 '16 at 22:36
  • 1
    I think you should explicitly put f[x,u[t,x]] in the latex formulation to make it clear if that's what you mean. – george2079 Dec 14 '16 at 22:39
  • @george2079 Thank you for your quick response, I edited the latex formulation. Now is it possible to solve it with a code that I wrote? – Hamed Peyrovedin Dec 14 '16 at 22:45
  • You should not delete your code. Please restore it. – Michael E2 Dec 08 '17 at 01:23

1 Answers1

4

The following illustrates how to use WhenEvent, but I had to change the NeumannValue at x == 0.05 to produce a well-behaved solution.

tmax = .01;
heat1 = NDSolveValue[{D[u[t, x], t] - (0.0000001 D[u[t, x], x, x] + b[t]*Exp[a[t] *x]) == 
    NeumannValue[0, x == 0] + NeumannValue[0 , x == 1/20],
    u[0, x] == 25, a[0] == 255 , b[0] == 2.5,
    WhenEvent[u[t, 1/20] > 250, {a[t] -> 0.01*a[t], b[t] -> .2 b[t]}]}, 
    {u, a, b}, {t, 0, tmax}, {x, 0, 1/20}, 
    Method -> {"MethodOfLines", "TemporalVariable" -> t, 
     "SpatialDiscretization" -> {"FiniteElement", 
    "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.001}}}},
    DiscreteVariables -> {a, b}];
Plot[Evaluate[heat1[[1]][t, 1/20]], {t, 0, tmax}, PlotRange -> All, 
    AxesLabel -> {t, "T(.05,t)"}]
LogPlot[heat1[[2]][t], {t, 0, tmax}, PlotRange -> All, AxesLabel -> {t, "a"}]

enter image description here

enter image description here

The second plot, added for convenience, shows the history of a.

The revised question has a rather different answer.

tmax = .001;
heat2 = NDSolveValue[{D[u[t, x], t] == 0.0000001 D[u[t, x], x, x] + 
    Piecewise[{{2.5 Exp[255 x], u[t, x] < 250}}, .5 Exp[2.55 x]] , 
    u[0, x] == 25, (D[u[t, x], x] /. x -> 0) == 0, 
    (D[u[t, x], x] /. x -> 1/20) == 7 (u[t, 1/20] - 25)}, 
     u, {t, 0, tmax}, {x, 0, 1/20}];
Plot[heat2[t, 1/20], {t, 0, tmax}, PlotRange -> All, AxesLabel -> {t, "T(.05,t)"}]

enter image description here

Perhaps more illustrative is

Plot3D[Evaluate[heat2[t, x]], {t, 0, tmax}, {x, 0, 1/20}, PlotRange -> All]

enter image description here

Note that the "FiniteElement" Method has been deleted, because it is incompatible as implemented in Mathematica 11.0.1 with the use of Piecewise, If, etc.

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156