André Schlichting

Classic results for porous medium equation

23 October 2012

Selfsimilar solutions

Let us seek for solutions of ∂ϱ=Δϱm\partial \varrho = \Delta \varrho^m satisfying the scaling hypothesis

(x,t)↦(y,τ):=(xs(t),τ(t))ϱ(x,t)=1s(t)Nu(xs(t),τ(t)), (x,t)\mapsto (y,\tau):=\left(\frac{x}{s(t)},\tau(t)\right) \qquad \varrho(x,t) = \frac{1}{s(t)^N} u\left(\frac{x}{s(t)},\tau(t)\right),

where s:R+→R+s: \R_+ \to \R_+ and τ:R+→R+\tau:\R_+ \to \R_+ are reparametrization of time and u:R+×RN→Ru:\R_+\times \R^N \to \R . The prefactor 1s(t)N\frac{1}{s(t)^N} ensures that uu conserves mass, i.e.

1=∫RNϱ(x,t)  dx=x↦ys∫RN1s(t)Nϱ(ys(t),t)  dy=∫RNu(y,τ(t))  dy. 1 = \int_{\R^N} \varrho(x,t) \; \dx{x} \stackrel{x\mapsto ys}{=} \int_{\R^N} \frac{1}{s(t)^N} \varrho(y s(t),t)\;\dx{y} = \int_{\R^N} u(y,\tau(t)) \; \dx{y}.

The time derivative has to satisfy

∂tϱ=τ′sN∂τu−Ns′sN+1u−s′sN+2x⋅∇yu=τ′s∂τu−s′sN+1∇y⋅(yu). \partial_t \varrho = \frac{\tau'}{s^N} \partial_\tau u - \frac{N s'}{s^{N+1}} u - \frac{s'}{s^{N+2}} x \cdot \nabla_y u = \frac{\tau'}{s} \partial_\tau u - \frac{s'}{s^{N+1}} \nabla_y \cdot ( y u) .

Further, let us calculate the gradient of ϱm\varrho^m

∇xϱ(x,t)m=∇y(1sNmu(y,τ)m)1s \nabla_x \varrho(x,t)^m = \nabla_y\left(\frac{1}{s^{Nm}} u(y,\tau)^m \right) \frac{1}{s}

and the Laplacian evaluates to

Δϱ(x,t)m=1sNm+2Δyu(y,τ)m. \Delta \varrho(x,t)^m = \frac{1}{s^{Nm+2}} \Delta_y u(y,\tau)^m .

Hence, the function uu solves the equation

τ′sN∂τu(y,τ)−s′sN+1∇y⋅(y  u(y,τ))=1sNm+2Δyu(y,τ)m. \frac{\tau'}{s^N} \partial_\tau u(y,\tau) - \frac{s'}{s^{N+1}} \nabla_y \cdot \bigl(y\; u(y,\tau)\bigr) = \frac{1}{s^{Nm+2}} \Delta_y u(y,\tau)^m .

We want the coefficients to be time-independent and a comparison results in the condition

c1τ′sN=c2s′sN+1=c31sNm+2=1, c_1 \frac{\tau'}{s^N} = c_2 \frac{s'}{s^{N+1}} = c_3 \frac{1}{s^{Nm+2}} = 1,

where c1,c2,c3>0c_1, c_2,c_3>0 are constants, which can be specified later.
From the first equality, we obtain τ′=c2c1s′s\tau' = \frac{c_2}{c_1}\frac{s'}{s} , which is solved by τ(t)=c2c1log⁡s(t)\tau(t) = \tfrac{c_2}{c_1}\log s(t) . The second equality, leads to

c3=s′c1sm(N−1)+1=c2ddts(t)N(m−1)+2N(m−1)+2 c_3 = s' c_1 s^{m(N-1)+1} = c_2 \frac{\dx}{\dx t} \frac{s(t)^{N(m-1)+2} }{N(m-1)+2}

and integrates to s(t)N(m−1)+2=c3c2(N(m−1)+2)t+c4s(t)^{N(m-1)+2} = \frac{c_3}{c_2} (N(m-1)+2)t+c_4 , where c4∈Rc_4\in \R is a further constant. Hence, we find the scaling relation

s(t)=(c3c2αt+c4)αandτ(t)=c2c1log⁡s(t),whereα:=1N(m−1)+2. s(t) = \left(\frac{c_3}{c_2\alpha} t+c_4\right)^{\alpha} \quad\text{and}\quad \tau(t) = \frac{c_2}{c_1} \log s(t) , \qquad\text{where}\qquad \alpha:= \frac{1}{N(m-1)+2} .

We are still free to choose the constants c1,c2,c3>0c_1,c_2,c_3>0 and c3∈Rc_3\in \R . A particular nice choice is given bys c1=1c_1=1 , c2=1αc_2=\frac{1}{\alpha} , c3=1c_3=1 and c4=0c_4=0 , then we obtain the result: If ϱ(x,t)\varrho(x,t) is a solution of the PME, then uu solves

∂τu(y,τ)=Δyum(y,τ)+α∇y⋅(y  u(y,τ)),withτ(t)=log⁡tands(t)=tα, \partial_\tau u(y,\tau) = \Delta_y u^m(y,\tau) + \alpha \nabla_y\cdot (y \; u(y,\tau)) , \quad\text{with}\quad \tau(t)=\log t \quad\text{and}\quad s(t) = t^\alpha ,

Equilibrium solutions

From the self similar rescaled solution uu , we can derive the equilibrium solution. Stationary solutions are given by function ϱ^:RN→R+\hat\varrho:\R^N\to \R_+ satisfying

∇y⋅(∇yϱ^m+αyϱ^)=0. \nabla_y \cdot\left(\nabla_y \hat\varrho^m + \alpha y \hat\varrho\right) = 0 .

Hence, by setting the flux inside of the divergence equal to zero

mϱ^m−1∇yϱ^+α  y  ϱ^=0. m \hat\varrho^{m-1} \nabla_y \hat\varrho + \alpha \;y\; \hat\varrho = 0 .

Hence, ϱ^≡0\hat\varrho\equiv 0 is a trivial solution and in the case m=1m=1 it is easy to check that ϱ^=e−α2y2\hat\varrho=e^{-\frac{\alpha}{2} y^2} is a solution (Compare this with the Ornstein-Uhlenbeck process, which is a special case of the Fokker-Planck equation}. Therefore, let us assume, that m≠1m\ne 1 . Then, we have

mϱ^m−2∇yϱ^+α  y=0, m \hat\varrho^{m-2} \nabla_y \hat\varrho +\alpha \;y = 0 ,

which can be rewritten as

∇y(mm−1ϱ^m−1)=∇y(−α2y2), \nabla_y \left(\frac{m}{m-1} \hat\varrho^{m-1}\right) = \nabla_y\left(- \frac{\alpha}{2} y^2\right),

which determines ϱ^m−1\hat\varrho^{m-1} up to a constant λ~∈R\tilde \lambda\in \R

ϱ^m−1=−m−1mα2y2+λ~ \hat\varrho^{m-1} = -\frac{m-1}{m} \frac{\alpha}{2}y^2 + \tilde \lambda

We can only take the power 1m−1\frac{1}{m-1} if the right hand side is non-zero, hence we set

ϱ^(y):={(λ~−m−1mα2y2)+1m−1, for m>1exp⁡(λ~−α2y2), for m=1(λ~−m−1mα2y2)1m−1, for m \hat\varrho(y) := \begin{cases} \left(\tilde \lambda- \frac{m-1}{m} \frac{\alpha}{2} y^2\right)_+^{\frac{1}{m-1}} &, \text{ for } m>1 \\ \exp\left(\tilde \lambda - \frac{\alpha}{2} y^2 \right) &, \text{ for } m = 1 \\ \left(\tilde \lambda - \frac{m-1}{m} \frac{\alpha}{2} y^2\right)^{\frac{1}{m-1}} &, \text{ for } m \end{cases}

hereby (x)+:=max⁡{0,x}(x)_+ := \max\{0, x\} denotes the positive part of xx . The constant λ~\tilde\lambda is chosen such that

∫RNϱ^(y)  dy=1. \int_{\R^N} \hat\varrho(y) \; \dx{y} = 1.