Problem with Piecewise Initial Conditions in NDSolve

Questions
, , ,
1

I’m trying to solve the system with piecewise initial conditions.

w = NDSolve[{D[u[r, thet, t], t] == 
    D[u[r, thet, t], {r, 2}] + (1/r) D[u[r, thet, t], r] + (1/r^2) D[
       u[r, thet, t], {thet, 2}] + 
     u[r, thet, t] (1 - u[r, thet, t] - c*x[r, thet, t]),
   D[x[r, thet, t], t] == 
    mu*(D[x[r, thet, t], {r, 2}] + (1/r) D[x[r, thet, t], 
          r] + (1/r^2) D[x[r, thet, t], {thet, 2}]) + 
     a*x[r, thet, t] (1 - c*u[r, thet, t] - b*x[r, thet, t]),
   D[u[r, thet, t], r] == 0 /. {r -> 1},
   D[x[r, thet, t], r] == 0 /. {r -> 1},
   
   D[u[r, thet, t], r] == 0 /. {r -> 0.01},
   D[x[r, thet, t], r] == 0 /. {r -> 0.01},
   (*u[r,thet,0]==0(*Piecewise[{{1,r<= r0},{0,r>r0}}]*),
   x[r,thet,0]==0(*Piecewise[{{1,r<= r0},{0,r>r0}}]*)},*)
   
   u[r, thet, t] == 
     Piecewise[{{Piecewise[{{0, 0 <= thet <= Pi/6}, {1, 
           Pi/6 < thet < 2 Pi}}], r <= r0}, {0, r > r0}}] /. {t -> 
      0},
   x[r, thet, t] == 
     Piecewise[{{Piecewise[{{1, 0 <= thet <= Pi/6}, {0, 
           Pi/6 < thet < 2 Pi}}], r <= r0}, {0, r > r0}}] /. {t -> 0}},
  {u[r, thet, t], x[r, thet, t]}, {t, 0, 5}, {r, 0.01, 1}, {thet, 0, 
   2 Pi}]

With this code Mathematica throws some errors

NDSolve::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable r.
NDSolve::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable thet. 
General::stop: "Further output of NDSolve will be suppressed during this calculation. "

How can I work this out? Without piecewise conditions Mathematica solves this system without problems.
Thanks.

2

Here are other ways which I tried to use:

Change boundaries of piecewise function like this

u[r, thet, 0] == 
  Piecewise[{{Piecewise[{{0, thet <= Pi/6}, {1, Pi/6 < thet}}], 
     r <= r0}, {0, r > r0}}],
 v[r, thet, 0] == 
  Piecewise[{{Piecewise[{{1, thet <= Pi/6}, {0, Pi/6 < thet}}], 
     r <= r0}, {0, r > r0}}]},

Or use NumericQ in piecewise function

foo1[r_?NumericQ, thet_?NumericQ] := 
 Piecewise[{{Piecewise[{{0, 0 <= thet <= Pi/6}, {1, 
       Pi/6 < thet < 2 Pi}}], r <= r0}, {0, r > r0}}]
foo2[r_?NumericQ, thet_?NumericQ] := 
 Piecewise[{{Piecewise[{{1, 0 <= thet <= Pi/6}, {0, 
       Pi/6 < thet < 2 Pi}}], r <= r0}, {0, r > r0}}]

The same error occurs

3

I tried to change initial function, like this, but get these errors again

b = 1; 
a = 1.01;
mu = 1;
r0 = 0.1;
r01 = 0.3;
r02 = 0.7;
r1 = 1;
foo1[r_?NumericQ] := 
 Piecewise[{{0, r < r0}, {(r - r0)/(r01 - r0), r0 <= r <= r01}, {1, 
    r01 < r < r02}, {(r - r1)/(r02 - r1), r02 <= r <= r1}, {0, 
    r > r1}}]


foo2[r_?NumericQ] := 
 Piecewise[{{1, r < r0}, {(r - r01)/(r0 - r01), r0 <= r <= r01}, {0, 
    r01 < r}}]


ftu = Table[{r, foo1[r]}, {r, 0, 1, 0.0025}];
ftv = ft = Table[{r, foo2[r]}, {r, 0, 1, 0.0025}];


initU = Interpolation[ftu];
initV = Interpolation[ftv];


w = NDSolve[{D[u[r, thet, t], t] == 
     D[u[r, thet, t], {r, 2}] + (1/r) D[u[r, thet, t], r] + (1/r^2) D[
        u[r, thet, t], {thet, 2}] + 
      u[r, thet, t] (1 - u[r, thet, t] - c*v[r, thet, t]),
    D[v[r, thet, t], t] == 
     mu*(D[v[r, thet, t], {r, 2}] + (1/r) D[v[r, thet, t], 
           r] + (1/r^2) D[v[r, thet, t], {thet, 2}]) + 
      a*v[r, thet, t] (1 - c*u[r, thet, t] - b*v[r, thet, t]),
    Derivative[1, 0, 0][u][1, thet, t] == 0, 
    Derivative[1, 0, 0][v][1, thet, t] == 0, 
    Derivative[1, 0, 0][u][0.01, thet, t] == 0, 
    Derivative[1, 0, 0][v][0.01, thet, t] == 0,
    (*D[u[r,thet,t],r]\[Equal]D[initU[r],r]/.{r\[Rule]1},
    D[v[r,thet,t],r]\[Equal]D[initV[r],r]/.{r\[Rule]1},
    D[u[r,thet,t],r]\[Equal]D[initU[r],r]/.{r\[Rule]0.01} ,
    D[v[r,thet,t],r]\[Equal]D[initV[r],r]/.{r\[Rule]0.01},*)
    
    u[r, 0, t] == u[r, 2 Pi, t],
    v[r, 0, t] == v[r, 2 Pi, t],
    u[r, thet, 0] == initU[r],
    v[r, thet, 0] == initV[r]},
   {u, v}, {t, 0, 1}, {thet, 0, 2 Pi}, {r, 0.01, 1}];

Can anybody give some advice?

4

Maybe there are other ways to solve this problem? I’ve got no more ideas. Can anybody help?

5

What happens if you restate your boundary conditions in the more conventional way of saying what the values are along the appropriate number of boundaries? Can you eliminate your nested Piecewise trying to specify your rectangle if you only attempt to solve over a rectangle of your space where your functions become constant for appropriate constant values of the parameters on the boundaries that make up the needed number of faces?

6

Thanks, Bill.
Do you mean something like this?

foo1[r_?NumericQ] := Piecewise[{{0, r < r0}, {1, r0 <= r <= 1}}, 0]
foo2[r_?NumericQ] := Piecewise[{{1, 0 <= r < r0}, {0, r >= r0}}]
ftu = Table[{r, foo1[r]}, {r, 0, 1, 0.0025}];
ftv = Table[{r, foo2[r]}, {r, 0, 1, 0.0025}];
initU = Interpolation[ftu];
initV = Interpolation[ftv];
w = NDSolve[{D[u[r, thet, t], t] == 
    D[u[r, thet, t], {r, 2}] + (1/r) D[u[r, thet, t], r] + (1/r^2) D[
       u[r, thet, t], {thet, 2}] + 
     u[r, thet, t] (1 - u[r, thet, t] - c*v[r, thet, t]),
   D[v[r, thet, t], t] == 
    mu*(D[v[r, thet, t], {r, 2}] + (1/r) D[v[r, thet, t], 
          r] + (1/r^2) D[v[r, thet, t], {thet, 2}]) + 
     a*v[r, thet, t] (1 - c*u[r, thet, t] - b*v[r, thet, t]),
   u[0.01, thet, t] == thet,
   u[1, thet, t] == thet,
   u[r, 0, t] == r,
   u[r, 2 Pi, t] == r,
   v[0.01, thet, t] == thet,
   v[1, thet, t] == thet,
   v[r, 0, t] == r,
   v[r, 2 Pi, t] == r,
   u[r, thet, 0] == foo1[r],
   v[r, thet, 0] == foo2[r]},
  {u, v}, {t, 0, 1}, {thet, 0, 2 Pi}, {r, 0.01, 1}, 
  PrecisionGoal -> 2]
7

Closer, but I was really hoping for something more like

r0 = 1/2;
foo1 = 0;
foo2 = 1;
w = NDSolve[{D[u[r,thet,t],t] == D[u[r,thet,t],{r,2}]+(1/r) D[u[r,thet,t],r]+(1/r^2) D[u[r,thet,t],{thet,2}] + 
    u[r,thet,t](1-u[r,thet,t]-c*v[r,thet,t]),
    D[v[r,thet,t],t] == mu*(D[v[r,thet,t],{r,2}]+(1/r)D[v[r,thet,t],r]+(1/r^2) D[v[r, thet, t], {thet, 2}]) +
    a*v[r, thet, t] (1 - c*u[r, thet, t] - b*v[r, thet, t]), 
   u[0.01, thet, t] == thet,
   u[1, thet, t] == thet,
   u[r, 0, t] == r,
   u[r, 2 Pi, t] == r,
   v[0.01, thet, t] == thet,
   v[1, thet, t] == thet,
   v[r, 0, t] == r,
   v[r, 2 Pi, t] == r,
   u[r, thet, 0] == foo1,
   v[r, thet, 0] == foo2},
  {u, v}, {t, 0, 1}, {thet, 0, 2 Pi}, {r, 0.01, r0}, PrecisionGoal -> 2]

which has hopefully removed anything confusing about the boundary conditions.

That results in

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>
NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

which strongly hints to me that there is a zero in a denominator somewhere, either explicitly or resulting from the processing of your derivatives or you are using a variable without having assigned it a constant value..

That leads me to discovering you have used mu without assigning a value to it. I assign 1 to mu and try again. That results in

NDSolve::bcedge: Boundary condition u[1,thet,t]==thet is not specified on a single edge of the boundary of the computational domain. >>

And I notice you have used a,b,c without assigning any values to those eiter. I assign 1 to each of those. With NDsolve you can’t use variables that haven’t been assigned values unless they are the functions you are solving for or the parameters for those.

That still gives me the same error. But this error is perhaps more specific and easier to track down. Can you make any progress on resolving this?

If I’d been able to solve the first problem then I would have followed up with

r0 = 1/2;
foo1 = 1;
foo2 = 0;
mu = 1;
a = 1;
b = 1;
c = 1;
w = NDSolve[{D[u[r,thet,t],t] == D[u[r,thet,t],{r,2}]+(1/r)D[u[r,thet,t],r]+(1/r^2)D[u[r,thet,t],{thet,2}] + 
     u[r,thet,t](1-u[r,thet,t]-c*v[r,thet,t]), 
   D[v[r,thet,t],t] == mu*(D[v[r,thet,t],{r,2}]+(1/r)D[v[r,thet,t],r]+(1/r^2)D[v[r,thet,t],{thet,2}]) + 
     a*v[r,thet,t](1-c*u[r,thet,t]-b*v[r,thet,t]), 
   u[0.01, thet, t] == thet,
   u[1, thet, t] == thet,
   u[r, 0, t] == r,
   u[r, 2 Pi, t] == r,
   v[0.01, thet, t] == thet,
   v[1, thet, t] == thet,
   v[r, 0, t] == r,
   v[r, 2 Pi, t] == r,
   u[r, thet, 0] == foo1,
   v[r, thet, 0] == foo2},
  {u, v}, {t, 0, 1}, {thet, 0, 2 Pi}, {r, r0, 1}, PrecisionGoal -> 2]

and then hoped I might have been able to stitch together the two solutions. Surprisingly, this actually gives a solution, but it still complains about inconsistent conditions. But this hints I may have gotten away from your zero. The warning is probably based on the same problem as the first case. If you can find and fix that you may be close.

Can you make any sense of this and make any more progress? Or perhaps someone else can see some problem with your conditions and point out how to fix this.

8

I have gotten this close

In[1]:= r0 = 1/2; foo1 = 0; foo2 = 1; mu = 1; a = 1; b = 1; c = 1;
w = NDSolve[{D[u[r, thet, t], t] == D[u[r, thet, t], {r, 2}] + (1/r) D[u[r, thet, t], r] + (1/r^2) D[
    u[r, thet, t], {thet, 2}] + u[r, thet, t] (1 - u[r, thet, t] - c*v[r, thet, t]), 
   D[v[r, thet, t], t] == mu*(D[v[r, thet, t], {r, 2}] + (1/r) D[v[r, thet, t], r] +
    (1/r^2) D[v[r, thet, t], {thet, 2}]) + a*v[r, thet, t] (1 - c*u[r, thet, t] - b*v[r, thet, t]),
   u[0.01, thet, t] == thet,
   u[r0, thet, t] == thet,
   u[r, 0, t] == r,
   u[r, 2 Pi, t] == r,
   v[0.01, thet, t] == thet,
   v[r0, thet, t] == thet,
   v[r, 0, t] == r,
   v[r, 2 Pi, t] == r,
   u[r, thet, 0] == foo1,
   v[r, thet, 0] == foo2}, {u, v}, {t, 0, 1}, {thet, 0, 2 Pi}, {r, 0.01, r0}, PrecisionGoal -> 2]

During evaluation of In[1]:= NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

Out[8]= {{u -> InterpolatingFunction[{{0.01, 0.5}, {0.,6.28318}, {0.,1.}}<>}, 
          v -> InterpolatingFunction[{{0.01, 0.5}, {0..6.28318}, {0.,1.}}<>]}}

In[9]:= r0 = 1/2; foo1 = 1; foo2 = 0; mu = 1; a = 1; b = 1; c = 1;
w = NDSolve[{D[u[r, thet, t], t] == D[u[r, thet, t], {r, 2}] + (1/r) D[u[r, thet, t], r] + (1/r^2) D[
    u[r, thet, t], {thet, 2}] + u[r, thet, t] (1 - u[r, thet, t] - c*v[r, thet, t]), 
   D[v[r, thet, t], t] == mu*(D[v[r, thet, t], {r, 2}] + (1/r) D[v[r, thet, t], r] +
    (1/r^2) D[v[r, thet, t], {thet, 2}]) + a*v[r, thet, t] (1 - c*u[r, thet, t] - b*v[r, thet, t]),
   u[r0, thet, t] == thet,
   u[1, thet, t] == thet,
   u[r, 0, t] == r,
   u[r, 2 Pi, t] == r,
   v[r0, thet, t] == thet,
   v[1, thet, t] == thet,
   v[r, 0, t] == r,
   v[r, 2 Pi, t] == r,
   u[r, thet, 0] == foo1,
   v[r, thet, 0] == foo2}, {u, v}, {t, 0, 1}, {thet, 0, 2 Pi}, {r, r0, 1}, PrecisionGoal -> 2]

During evaluation of In[9]:= NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

Out[16]= {{u -> InterpolatingFunction[{{0.5, 1.},{0.,6.28318},{0., 1.}}<>]
           v -> InterpolatingFunction[{{0.5, 1.},{0.,6.28318},{0., 1.}}<>]}}

Is it possible that two of your conditions do not match at one of the corners? Yes. Take every pair of your conditions and see if you can force the two left hand sides to be the same and the right hand sides to be different. When you have the list of all the inconsistencies then see if you can repair these. Then we see if the two solutions can be made compatible along the boundary r0.

9

Thank you, Bill.