# Solve PDE with complicated coefficient non-linearity [closed]

I wish to solve is the heat equation with solution-dependent coefficient.
The equation along with BC's and IC's are as follows:
$$\begin{equation} \frac{\partial P[x,t]}{\partial t}-\alpha[x,t]\frac{\partial^2 P[x,t]}{\partial x^2}=0 \end{equation}$$,
$$\begin{equation}P[x,0]=0,P[L,t]=0 \end{equation}$$,
$$P[0,t]=Step[t-t1]-0.5Step[t-t2]-0.25Step[t-t3]$$,
$$L=1$$.
The $$\alpha[x,t]$$ coeffiecnt is determined by the following rules:
$$\alpha[x,0]=\alpha_1$$,
If $$P[x,t]>0.7$$ $$\rightarrow$$ Set $$\alpha[x,t]$$ to $$\alpha_2$$
If $$P[x,t]<0.3$$ $$\rightarrow$$ Set $$\alpha[x,t]$$ to $$\alpha_1$$
If $$0.3 $$\rightarrow$$ Do Nothing.

Those conditions represent the case of bistability of the diffusion coefficient.

I've tried quite a few methods to solve it, including the definition of $$\alpha$$ as a function of $$\frac{\partial P[x,t]}{\partial t}$$ and $$P[x,t]$$ with a combination of Sing[],Boole[],WhenEvent[],ect.

None of them worked at the past few weeks, so any idea may help.
Thanks,
Ofek.

• If the thermal conductivity is variable, don't you need to express the heat equation as $\rho {{\hat C}_p}(u)\frac{{\partial u}}{{\partial t}} + \nabla \cdot \left( { - k(u)\nabla u} \right) = 0$? Jan 5 '19 at 13:49
• Yes, the case I've mentioned is a simplified version of the full equation (as you wrote) Jan 5 '19 at 14:14
• In the last application of the WhenEvent Documentation, they have an example of an ON/OFF controller (which is a bistable system with hysteresis) applied to the 2D heat equation. Did you try to adapt that example? Jan 5 '19 at 15:59
• @OfekPeretz Can you provide data for a numerical model? Jan 5 '19 at 16:26
• It is effectively impossible to diagnose what has gone wrong if the actual code used is not posted. Jan 5 '19 at 16:26

I think the OP is trying to demonstrate hysteresis. WhenEvent has an example of an on-off controller, which is an example of a bistable system. Unfortunately, my understanding is, that an event must be a function of time only and cannot have local spatial dependence. The essence of the Numerical Method of Lines is to turn a time-dependent PDE into a system of ODEs depending on time only. This suggests that perhaps we can combine WhenEvent with a roll your own version of the Numerical Method Of Lines. To build this model, it is best to start simple so let us start with a zero dimensional case that demonstrates hysteresis with WhenEvent.

## 0-D Case

Let us consider a small system with no diffusion limitations driven by a fluctuating heat source. The ODE that we wish to solve is given by:

$${\hat C_p}(u)\frac{du}{dt} = q\left( t \right)$$

We will start with a system that is initially solid. We will set up a simple test case to verify that our WhenEvent construction gives us the correct bistable behavior.

cpl = 2;
cps = 1;
thresholdliquid = Rationalize[0.7];
thresholdsolid = Rationalize[0.3];
q[t] = 5/2 UnitStep[-7 + t] - 5/2 UnitStep[-6 + t] +
5/2 UnitStep[-5 + t] - 5/2 UnitStep[-3 + t] + UnitStep[-1 + t];
NDSolve[{cp[t] u'[t] == q[t], cp == cps, u == 0,
WhenEvent[{u[t] == thresholdliquid,
u[t] == thresholdsolid} , {cp[t] ->
If[u[t] == thresholdliquid, cpl,
If[u[t] == thresholdsolid, cps, cp[t]]]}]}, {u, cp}, {t, 0,
10}, DiscreteVariables -> {cp \[Element] {cps, cpl}}  ];
Plot[{Evaluate[{u[t], cp[t] - cps} /. %], thresholdsolid,
thresholdliquid}, {t, 0, 10}, PlotTheme -> "Web",
PlotStyle -> {Thick, Thick, Dashed, Dashed}] You can clearly see breaks in the temperature curve that have different thresholds if the temperature is rising or falling. So, it looks as though we have a good prototype for bistability using WhenEvent in the 0D case.

## Combining WhenEvent with a Roll-Your-Own Numerical Method of Lines

In the Boundary Conditions section of The Numerical Method Of Lines Tutorial, Shows how you can roll your own numerical method of lines schemes to convert a PDE into a system of ODEs. The following Mathematica code attempts to combine our WhenEvent prototype with the procedure outlined in the numerical method of lines tutorial. I used arbitrary values for the parameters since none were supplied. The thermal diffusivity difference between ice and water is about a factor of 9. However, an ice water system would have such a large latent heat effects that you would need to use the form of the equation I alluded to in the comments.

(* Set Parameters *)
n = 50; Subscript[h, n] =
1/n; sf = 1; tmax = 10; alphasolid = 1; alphaliquid = 1/10*alphasolid;
t1 = 1; t2 = 4; t3 = 5 ;
thresholdliquid = Rationalize[0.7];
thresholdsolid = Rationalize[0.3];
(* Create Temperature and Phase State Vector *)
U[t_] = Table[Subscript[u, i][t], {i, 0, n}];
A[t_] = Table[Subscript[a, i][t], {i, 0, n}];
dvbls = Head[#] & /@ A[t];
(* Differentiate BCs *)
bc0 = Subscript[u, 0][t] ==
UnitStep[t - t1] - 1/2*UnitStep[t - t2] - 1/4*UnitStep[t - t3];
bc0d = Map[D[#, t] + sf # &, bc0];
bc1 = Subscript[u, n][t] == 0;
bc1d = Map[D[#, t] + sf # &, bc1];
(* Set Laplacian *)
laplacianu =
NDSolveFiniteDifferenceDerivative[2, Subscript[h, n] Range[0, n],
U[t]];
(* Set up system of ODEs *)
eqns = Thread[D[U[t], t] -  A[t] laplacianu == 0];
(* Replace Boundary Condtions Equations *)
eqns[] = bc0d;
eqns[[-1]] = bc1d;
(* Set Initial Temperature to 0 *)
ics = Thread[U == Table[0, {n + 1}]];
(* Set Initial Phase State to Solid *)
istates = Thread[A == Table[alphasolid, {n + 1}]];
(* Set up Events*)
events = WhenEvent[{Subscript[u, #][t] == tl,
Subscript[u, #][t] == ts} , {Subscript[a, #][t] ->
If[Subscript[u, #][t] == tl, alphaliquid,
If[Subscript[u, #][t] == ts, alphasolid,
Subscript[a, #][t]]]}] & /@
Range[0, n] /. {ts -> thresholdsolid, tl -> thresholdliquid};
(* Join System of Equations to be Solved *)
system = Join[eqns, ics, istates, events];
lines = NDSolve[system, Join[U[t], A[t]], {t, 0, tmax},
DiscreteVariables -> dvbls];
(* Create Interpolation Functions for Easier Plotting *)
f = Interpolation[
Flatten[Table[{{t, i Subscript[h, n]},
First[Subscript[u, i][t] /. lines]}, {i, 0, n}, {t, 0, tmax,
tmax/500}], 1], InterpolationOrder -> 2];
fa = Interpolation[
Flatten[Table[{{t, i Subscript[h, n]},
First[Subscript[a, i][t] /. lines]}, {i, 0, n}, {t, 0, tmax,
tmax/500}], 1], InterpolationOrder -> 2];
(* Plot Temperature and Diffusivity *)
ContourPlot[f[t, x], {t, 0., tmax}, {x, 0., 1},
PlotPoints -> {200, 200}, MaxRecursion -> 4, PlotTheme -> "Web",
ColorFunction -> ColorData["DarkBands"]]
ContourPlot[fa[t, x], {t, 0., tmax}, {x, 0., 1},
PlotPoints -> {200, 200}, MaxRecursion -> 4, PlotTheme -> "Web",
ColorFunction -> ColorData["DarkBands"]]


### Temperature Contour Plot ### Diffusivity State Contour Plot ### Demonstration of Bistability

We should be able to demonstrate bistability by plotting the temperature and the diffusivity at a value that intersects the phase envelope (e.g., $$x=0.1$$) as shown with the following Mathematica code.

Plot[{f[t, 0.1], fa[t, 0.1], thresholdsolid, thresholdliquid}, {t, 0,
10}, PlotTheme -> "Web",
PlotStyle -> {Thick, Thick, Dashed, Dashed}] The system behaves as expected showing transitions at 0.7 going up and 0.3 going down. One can use ParametricPlot to demonstrate the path dependence more clearly as shown below.

p0 = ParametricPlot[{f[t, 0.1], fa[t, 0.1]}, {t, 0, 10},
ColorFunction -> (ColorData["Rainbow"][#3] &), MaxRecursion -> 0,
PlotPoints -> 100, PlotTheme -> "Web"];
p0 /. l_Line :> {Arrowheads[ConstantArray[0.05, 30]], Arrow[l]} ## Update To Include No Flux Condition at x=1

In the comments, Ofek requested to also consider the case of a no flux condition at $$x$$=$$1$$. It turns out that the MOL Tutorial actually considers this very case. One just needs to adjust bc1 and bc1d accordingly

bc1 = Last[
NDSolveFiniteDifferenceDerivative[1, Subscript[h, n] Range[0, n],
U[t]]] == 0;
bc1d = Map[D[#, t] &, bc1];


Note that the scale factor $$sf$$ has been set to zero for Neumann type conditions. Without adjustments to the first case, you will obtain plots for temperature and diffusivity state as shown below:

### Temperature Contour Plot No Flux on End ### Diffusivity State Contour Plot No Flux on End You will have to adjust parameters to make the plots look more interesting.

• a very good example in MMA! Jan 7 '19 at 7:13
• @TimLaska, Thanks for your detailed solution, this is exactly what I wanted to solve. Just one more question - How can I change the boundary condition at $X=1$ from $P_n(0)==0$ into $\frac{\partial P_n(0)}{\partial x}==0$? Jan 14 '19 at 12:29
• @OfekPeretz I added a no flux case to the post. Jan 14 '19 at 13:46
• @TimLaska, Amazing. Thanks! Jan 14 '19 at 15:12
• @OfekPeretz, You are welcome! Jan 14 '19 at 15:23

I do not see what the problem is.

t1 = 1; t2 = 2; t3 = 3; L = 10; T = 10;
bc = {P[t, 0] == 0,
P[t, L] ==
UnitStep[t - t1] - UnitStep[t - t2]/2 - UnitStep[t - t3]/4};
a[p_] := Piecewise[{{1, p < .3}, {2, .3 <= p < .5}, {1/2, True}}]
ic = {P[0, x] == 0};
eq = D[P[t, x], t] - a[P[t, x]]*D[P[t, x], x, x] == 0;
sol = NDSolveValue[{eq, ic, bc}, P, {t, 0, T}, {x, 0, L}];
{ContourPlot[sol[t, x], {t, 0, T}, {x, 0, L}, Contours -> 20,
ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All,
FrameLabel -> {"t", "x"}],
ContourPlot[sol[t, x], {t, 0, 5}, {x, L - 2, L}, Contours -> 20,
ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All,
FrameLabel -> {"t", "x"}]} 