MATH529 L06
Heat equation

Relation to Brownian motion

Consider (Subscript[N, p]) initially within (Subscript[r, x])=[-1/2,1/2] executing random motion in 1D.

Np=10;rx={-1,1}/2; x=RandomReal[rx,Np]; ry=rx/50; r={rx,ry};
Graphics[Point[{#,0}]& /@ x,PlotRange->r,Axes->{True,False}]

Consider (Subscript[N, p])≫1 indistinguishable particles

atoms[n_]:=Graphics[Point[{#,0}]& /@ RandomReal[rx,n],PlotRange->r,Axes->{True,False}];

Manipulate[atoms[n],{n,100,1000,100}]

Density function

Rather than keeping track of individual particle, introduce a density function u(t,x; h) such number of particles in interval [x-h/2,x+h/2] is n(t,x; h)=h u(t,x; h)

Np=(10^6);x=RandomReal[rx,Np]; Nb=100; h=1/Nb; Histogram[x,Nb]

Particle density  evolution in time

Normal random particle velocity of mean μ=0 and variance (σ^2)

μ=0; σ=1; v[n_]:=RandomVariate[NormalDistribution[μ,σ],n];
Histogram[v[Np],Nb]

xt[Δt_] := x + Δt v[Np];

Manipulate[Histogram[xt[Δt],{h},PlotRange->{{-5,5},{0,(10^4)}}],{Δt,0,2,0.02}]

Particle density  evolution in time

Normal random particle velocity of mean μ=0 and variance (σ^2)

v[n_]:=RandomVariate[UniformDistribution[{-1/2,1/2}],n];
Histogram[v[Np],Nb]

xt[Δt_] := x + Δt v[Np];

Manipulate[Histogram[xt[Δt],{h},PlotRange->{{-5,5},{0,(10^4)}}],{Δt,0,2,0.02}]

Separation of variables

Try u(x,y)=X(x)Y(y) in (Subscript[u, x x])=4 (Subscript[u, y])⇔(Subscript[u, x x])- 4 (Subscript[u, y])=0, homogeneous PDE

PDE = (∂_(x,x))u[x,y]==4 (∂_y)u[x,y]; vPDE = (∂_(x,x))v[x,y]==4 (∂_y)v[x,y];
u[x_,y_]=X[x] Y[y];
sPDE=Assuming[u[x,y]≠0,DivideSides[PDE,u[x,y]]]

xDE = sPDE[[1]] == -λ; yDE = sPDE[[2]] == -λ; {xDE,yDE}

ODE solutions depend on sign of λ, consider each possible case

Case   λ=0

xDE0 = xDE /. λ->0; yDE0 = yDE /. λ->0;
{xDE0,yDE0}

xsol0=DSolve[xDE0,X[x],x][[1,1]]; ysol0=DSolve[yDE0,Y[y],y][[1,1]] /. (C[1])->(C[3]) ;
sol0={xsol0,ysol0}

u0[x_,y_] = X[x] Y[y] /. sol0

Verify solution

v[x_,y_] = u0[x,y];
vPDE

Clear[v]; vPDE

Case   λ<0

xDEn = xDE /. λ->-(α^2); yDEn = yDE /. λ->-(α^2);
{xDEn,yDEn}

xsoln=DSolve[xDEn,X[x],x][[1,1]]/. {(C[1])->(C[4]),(C[2])->(C[5])}; ysoln=DSolve[yDEn,Y[y],y][[1,1]] /.(C[1])->(C[6]) ;
soln={xsoln,ysoln}

un[x_,y_] = X[x] Y[y] /. soln

Verify solution

v[x_,y_] = un[x,y];
Simplify[vPDE ]

Clear[v]; vPDE

Case   λ>0

xDEp = xDE /. λ->(α^2); yDEp = yDE /. λ->(α^2);
{xDEp,yDEp}

xsolp=DSolve[xDEp,X[x],x][[1,1]]/. {(C[1])->(C[7]),(C[2])->(C[8])}; ysolp=DSolve[yDEp,Y[y],y][[1,1]] /.(C[1])->(C[9]) ;
solp={xsolp,ysolp}

up[x_,y_] = X[x] Y[y] /. solp

Verify solution

v[x_,y_] = up[x,y];
Simplify[vPDE ]

Clear[v]; vPDE

Superposition of solutions

v[x_,y_]= a0 u0[x,y] + ap up[x,y] + an un[x,y]

Simplify[vPDE ]