# Piecewise causes ndnum warning in NDSolve

Bug introduced in 8.0.4 or earlier.

I try to solve the heat transfer equation with boundary condition that depends on time:

r0 = 0.75 10^-3;(*Beam spot size, m*)
ω = π ν;
ν = 1; (*pulse repetition rate, Hz*)

c = 1710;(*Heat capacity, W/(m·K)*)
ρ = 879; (*Density, kg/m^3*)
λ = 0.111;(*Heat conductivity, W/(m·K)*)
T0 = 300;(*initial temperature,K*)
T1 = 1000 - T0; (*Hot state temperature, K*)
Rm = 3 10^-3;(*Sample Radius, m*)
zm = 2 10^-3 ;(*Sample thickness, m*)

eq = D[T[R, z, t], t] == λ/(
c ρ) (1/(R + 10^-20) D[R D[T[R, z, t], R], R] +
D[T[R, z, t], {z, 2}]);

init1 = T[R, z, 0] == T0;
bc1 = D[T[R, z, t], {R, 1}] == 0 /. R -> 0;
bc2 = D[T[R, z, t], {R, 1}] == 0 /. R -> Rm;
bc3 = T[R, z, t] ==
Piecewise[{{T0 + T1 (1 - Abs@Sin[ω t])^500,
0 <= R <= r0}, {T0, True}}] /. z -> 0;
bc4 = T[R, z, t] == T0 /. z -> zm;

sol = NDSolveValue[{eq, init1, bc1, bc2, bc3, bc4},
T[R, z, t], {R, 0, Rm}, {z, 0, zm}, {t, 0, 10},
AccuracyGoal -> 30, MaxStepFraction -> 0.05]


The piecewise function defines the temperature at the left side of the z-interval as a short square pulses coming with frequency ω to central part of the disk.

However, solver returns the error-messages:

NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0.

What's wrong? I've tried with FEM package but it produce even more error-messages :)

• When t=0, init and bc3 are inconsistent, it makes T[R, 0, 0] both 300 and 1000 for 0 <= R <= 0.00075. – Feyre Feb 15 '17 at 9:43
• @Feyre, ok, let's add small shift to t at bc3 making (t-5*10^-3) instead of t inside Sin. Unfortunately, it does not change principally anything. – Rom38 Feb 15 '17 at 10:42
• Why don't change Sin to Cos in bc3? – zhk Feb 15 '17 at 10:51
• @MMM, Is it principal? :) – Rom38 Feb 15 '17 at 10:53
• @Rom38 I don't know about that. – zhk Feb 15 '17 at 10:54

Several issues here.

1. As mentioned by Feyre, the i.c. and b.c. are inconsistent. I'll set "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100} in response. For more information check this post.

2. ndnum is caused by a bug of NDSolve (yeah it's confirmed, I reported it WRI before), similar problem has been observed in this and this post. This can be resolved by adding SimplifyPWToUnitStep@ before Piecewise.

3. The modeling for laser pulse is improper. Ideally DiracDelta is a possible choice, but currently NDSolve can't handle it properly, so we need to use an approximate one, for example:

dirac[r_, a_] = Sqrt[a/Pi] Exp[-a r^2];


Periodity can be achieved with Mod:

Plot[dirac[Mod[t, 1, -1/2], 100], {t, 0, 10}, PlotRange -> All]


Here's the solution:

(* Previous code is not modified so I'd like to omit them in this post. *)
dirac[r_, a_] = Sqrt[a/Pi] Exp[-a r^2];

bc3 = T[R, z, t] == SimplifyPWToUnitStep@
Piecewise[{{T0 + T1*dirac[Mod[t, Pi/ω, -2^(-1)], 100],
0 <= R <= r0}, {T0, True}}] /. z -> 0;
bc4 = T[R, z, t] == T0 /. z -> zm;

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf : False | True, sf_: Automatic] := {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}

sol = NDSolveValue[{eq, init1, bc1, bc2, bc3, bc4},
T, {R, 0, Rm}, {z, 0, zm}, {t, 0, 10},
Method -> Union[mol[60, 2], mol[True, 100]]]; // AbsoluteTiming

Animate[Plot3D[sol[r, z, t], {r, 0, Rm}, {z, 0, zm}, PlotRange -> {0, 4000},
PerformanceGoal -> "Quality"], {t, 0, 10}]


Remark

1. NDSolve spits out eerr warning, but I think it's not a big problem since the error is small. Use more grid points will probably suppress the error, you can have a try if you have time.

2. I think "FiniteElement" will do a better job on solving this problem.

3. I took away AccuracyGoal because of the reason mentioned here.

• Thanks a lot! I'm enough new in solution of DE in Mathematica.. – Rom38 Feb 16 '17 at 5:46