WhenEvent and partial derivatives

Question

Can WhenEvent be used to reset the conditions on a PDE at a given time? How would the syntax of that be?

This is the code I`m using

r = 2;
d2 = 1;
l = 5;
sol = NDSolve[{Derivative[0, 1][n][x, t] ==
               d2 Derivative[2, 0][n][x, t] + r*n[x, t]*(1 - n[x, t]),
               n[x, 0] == 0.1*x*(1 - x/l), n[0, t] == 0, n[l, t] == 0,
               WhenEvent[t == 50, n[x, t] -> 0]}, {n}, {x, 0, l}, {t, 0, tmax}]

Could you maybe specify the particular PDE you're working with, so that it's easier for us to show you how? — J. M.'s missing motivation, Jun 09 '13 at 04:29
Im trying to solve a reaction diffusion equation with NDSolve, and reset the density function mid integration — fernando, Jun 09 '13 at 04:36
sol = NDSolve[{ !( *SubscriptBox[([PartialD]), (t)](n[x, t])) == d2*!( *SubscriptBox[([PartialD]), (x, x)](n[x, t])) + rn[x, t](1 - n[x, t]), n[x, 0] == 0.1x(1 - x/l), n[0, t] == 0, n[l, t] == 0, WhenEvent[t == 50, n[x, t] -> 0]}, {n}, {x, 0, l}, {t, 0, tmax}, — fernando, Jun 09 '13 at 04:39
Well, please edit that code into your post, instead of posting that in comments... — J. M.'s missing motivation, Jun 09 '13 at 04:43
No problem, but without some typical values for the parameters (l, d2, etc.), it's hard to test. Could you add those to the question? — Michael E2, Jun 09 '13 at 05:08

Silvia · Answer 1 · 2013-06-10T18:58:22.730

17

According to the error message:

NDSolve::nlnum1: "The function value {0} is not a list of numbers with dimensions {25} when the arguments are {50.,{<<25>>}."

I think you should feed a 25-length list of 0 to n[x,t] in the WhenEvent:

WhenEvent[t > 50, n[x, t] -> ConstantArray[0, 25]]

Plot the result:

Plot3D[Evaluate[n[x, t] /. sol], {x, 0, l}, {t, 0, 100}, PlotPoints -> 50]

enter image description here

Edit:

According to the documentation, NDSolve automatically does processing for discontinuous functions like Sign, so here an alternative way which do not require manual specifying the number of grid nodes (i.e. 25):

sol2 = NDSolve[{
   Derivative[0, 1][n][x, t] ==
             d2*Derivative[2, 0][n][x, t] + 
             r*Sign[50 - t]*n[x, t]*(1 - Sign[50 - t]*n[x, t]),
    n[x, 0] == 0.1*x*(1 - x/l),
    n[0, t] == 0, 
       n[l, t] == 0},
   {n}, {x, 0, l}, {t, 0, 100}]

Note the difference near the discontinuity line between this solution and sol by above WhenEvent:

Show[
 MapThread[Plot3D[
    Evaluate[n[x, t] /. #1], {x, 0, l}, {t, 49, 53},
    PlotPoints -> 50, PlotStyle -> None,
    MeshFunctions -> (#1 &), MeshStyle -> #2,
    BoundaryStyle -> None, PlotRange -> All
    ] &,
  {{sol, sol2}, {Red, Blue}}
  ]]

enter image description here

Edit 2:

Using WhenEvent with automatic detecting the x-grid:

<< DifferentialEquations`InterpolatingFunctionAnatomy`
Clear[xGridExtractor]
xGridExtractor[f_] := InterpolatingFunctionCoordinates[Head[f]][[1]]

sol = NDSolve[
        {
         Derivative[0, 1][n][x, t] == d2*Derivative[2, 0][n][x, t] + r*n[x, t]*(1-n[x, t]), 
         n[x, 0] == 0.1*x*(1 - x/l),
         n[0, t] == 0,
         n[l, t] == 0, 
         WhenEvent[t > 50, n[x, t] -> 0*xGridExtractor[n[x, t]]]
        },
        {n}, {x, 0, l}, {t, 0, 100}]

Edit 3:

According to OP's comment, here is how to reset the initial condition along $t=50$ in a more general sense:

sol = NDSolve[
        {
         Derivative[0, 1][n][x, t] == d2*Derivative[2, 0][n][x, t] + r*n[x, t]*(1-n[x, t]), 
         n[x, 0] == 0.1*x*(1 - x/l),
         n[0, t] == 0,
         n[l, t] == 0, 
         WhenEvent[t > 50, n[x, t] -> (.5 Head[n[x, t]] /@ xGridExtractor[n[x, t]])]
        },
        {n}, {x, 0, l}, {t, 0, 100}]

Plot3D[Evaluate[n[x, t] /. sol], {x, 0, l}, {t, 49, 51},
     PlotPoints -> 50, MeshFunctions -> {#2 &, #3 &}, Exclusions -> t == 50]

half value

edited Jun 10 '13 at 18:58

answered Jun 09 '13 at 05:38

Silvia

27,556
3
84
164

+1. I was just writing such an answer! The array seems to reinitialize the values of n[x, 50] along the line 0 <= x <= l -- probably the interpolation grid. – Michael E2 Jun 09 '13 at 05:40
@MichaelE2 Thanks and agree. I think there should be an alternative way which does not require manual specifying the number of grid nodes. – Silvia Jun 09 '13 at 05:53
I like the use of Sign[], but maybe the PDE will look nicer if you use Piecewise[]. – J. M.'s missing motivation Jun 09 '13 at 20:26
@0x4A4D I just tried Piecewise, and I think I'd prefer the WhenEvent version. Neither Sign nor Piecewise are as good as WhenEvent near the discontinuity interface. – Silvia Jun 09 '13 at 20:36
That's too bad. I was hoping for an alternative that'd also work in earlier versions of Mathematica. (I'm still on version 8.) – J. M.'s missing motivation Jun 09 '13 at 20:41
@0x4A4D In that case I guess we'd have to increase the interpolation points on x-grid.. – Silvia Jun 09 '13 at 20:44
@MichaelE2 I've spent some time trying to adapt your solutions for a case in which the density is halfed instead of zeroed. Do you have any suggestions? – fernando Jun 10 '13 at 18:26
@MichaelE2 Just for my curiosity, how did you manage to @ more than one person? – Silvia Jun 10 '13 at 18:35
@MichaelE2, sorry for mistagging! – fernando Jun 10 '13 at 18:35
@Silvia I think he can do that because im the thread owner, so im automatically notified – fernando Jun 10 '13 at 18:36
@fernando Thanks to inform :) Do you mean you want n[x, t] -> 0.5 n[x, t]? – Silvia Jun 10 '13 at 18:39
@Silvia, yes, that's what I'm trying now! By the way, I feel like I owe you a cupcake or something for all the help, thanks a lot! – fernando Jun 10 '13 at 18:42
@fernando You are welcome :) Please see my edit 3. – Silvia Jun 10 '13 at 18:58
@MichaelE2 Thanks for informing. I thought there are some tricks :D – Silvia Jun 10 '13 at 19:20

WhenEvent and partial derivatives

1 Answers1

Linked