# NDSolve for a system of PDEs with piecewise coefficients

I want to solve a system of 3 second order linear PDEs with homogeneous Dirichlet boundary conditions. The code I use is as follows:

ν = 0.3;
ϰ = 0.5;
lh = 5;
dh = 0.1;
Na = 50;
Eib = 5;
Feb = 10^-6;

VB = lh*lh*1;
Vbd = 1/6 π dh^3;
VBd = VB - Na Vbd;
x01[n_] := 2 n dh
x02[n_] := 2 n dh
x03 = 0.5;
x1F = 2;
x2F = 2;

chib[x1_, x2_, x3_] := 1 - chii[x1, x2, x3]
chii[x1_, x2_, x3_] := Sum[UnitStep[0.25 dh^2 - (x1 - x01[n])^2 - (x2 - x02[m])^2 - (x3 - x03)^2], {n,1,24}, {m,1,24}]
chiF[x1_, x2_] := UnitStep[1^2 - (x1 - x1F)^2 - (x2 - x2F)^2]

F[x1_, x2_] := -Feb chiF[x1, x2]
Ye[x1_, x2_, x3_] :=
VB/VBd chib[x1, x2, x3] + Eib VB/Vbd chii[x1, x2, x3]
A0[x1_, x2_] := 1/(1 - ν^2) Integrate[Ye[x1, x2, x3], {x3, 0, 1}]
A2[x1_, x2_] :=
1/(1 - ν^2) Integrate[x3^2 Ye[x1, x2, x3], {x3, 0, 1}]

sol = NDSolve[{D[A2[x1, x2] (D[ϕ1[x1, x2], x1] + ν D[ϕ2[x1, x2], x2]), x1] + 0.5*(1 - ν) D[A2[x1, x2] (D[ϕ1[x1, x2], x2] + D[ϕ2[x1, x2], x1]), x2] + 0.5*ϰ (1 - ν) A0[x1, x2] (D[u3[x1, x2], x1] - ϕ1[x1, x2]) == 0, 0.5*(1 - ν) D[A2[x1, x2] (D[ϕ1[x1, x2], x2] + D[ϕ2[x1, x2], x1]), x1] + D[A2[x1, x2] (ν D[ϕ1[x1, x2], x1] + D[ϕ2[x1, x2], x2]), x2] + 0.5*ϰ (1 - ν) A0[x1, x2] (D[u3[x1, x2], x2] - ϕ2[x1, x2]) == 0, 0.5*ϰ (1 - ν) D[A0[x1, x2] (D[u3[x1, x2], x1] - ϕ1[x1, x2]), x1] + 0.5*ϰ (1 - ν) D[A0[x1, x2] (D[u3[x1, x2], x2] - ϕ2[x1, x2]), x2] + F[x1, x2] == 0, u3[0, x2] == 0, u3[lh, x2] == 0, u3[x1, 0] == 0, u3[x1, lh] == 0, ϕ1[0, x2] == 0, ϕ1[lh, x2] == 0, ϕ1[x1, 0] == 0, ϕ1[x1, lh] == 0, ϕ2[0, x2] == 0, ϕ2[lh, x2] == 0, ϕ2[x1, 0] == 0, ϕ2[x1, lh] == 0}, {u3, ϕ1, ϕ2}, {x1, 0, lh}, {x2, 0, lh}]


The programm runs for several hours and does not compute the solution. I suspect that the difficulty is in the definition of chii: when there is only 1 term, NDSolve gives the solution within several seconds. Probably, many-term UnitStep complicates the computations.

I have tried to manipulate with WorkingPrecision and with method. In particular, I have used

Method -> {"DiscontinuityProcessing" -> False}

Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}}

Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 1000}}}

Method -> {"DAEInitialization" -> {"Collocation", "CollocationDirection" -> "Forward"}}

StartingStepSize -> 0.1


but it computes several hours and does not give anything.

Any hint or suggestion is highly appreciated.

• It's not NDSolve that makes your code slow. Execute a single A0[x1, x2] and you'll know what I mean. Use NIntegrate and ?NumericQ to redefine A0 and A2 may help. (Haven't tested it, now I only have a old laptop at hand. ) Commented Aug 21, 2018 at 2:12
• When there is only one term in the definition of chii, NDSolve computes the solution in a minute. ?NumericQ and NIntegrate does not help much. Commented Aug 21, 2018 at 2:40
• That's simply because Integrate can compute the integral fast when chii isn't complicated. You can move the equations out of NDSolve and timing them separately. Commented Aug 21, 2018 at 2:47
• Maybe there is a way to redefine the characteristic function chii` to avoid heavy computations. Commented Aug 21, 2018 at 2:50