Numerical approximation of heat equation in bounded domains in Mathematica

Oana Rachieru

Transilvania University of Brasov

We consider numerical approximation of heat equation with Neumann condition and Dirichlet condition on the boundary of domain  D⊂(R^d) , that is we will consider the problem of finding a (C^(1, 2)) 

              function\(u = u(t, z) : \((0, ∞)\)×D\&_\)->R satisfying the Neumann condition on the  boundary domain, respectively Dirichlet condition on the boundary domain  .  
            Using probabilistic methods we can give explicit representation of solutions of heat equation on both problems.

Introduction of Brownian Motion

Brownian motion →observed by Robert Brown 1828

→iregular motion of small pollen grains on the surface of  a liquid
1923  -  
        N .  Wiener proved proved the consistency of the definiton with the specified properties
1951  -   D . M . 
            Donsker proved convergence of random walk to BM

Construction

BM  →  constructed as the limit of a random walk

consider space (Ω, F, P)

consider i . id random variables (Subscript[Y, n]) with:

\(P \((Y\_n = 1)\) = \(P \((Y\_n = \(-1\))\) = 1\/2\), \ 
      n ∈ \(ℕ\^*\)\)

and define the random walk (Subscript[S, n]) by

(Subscript[S, n])={({{0,  n = 0}, {(∑_(i = 1)^n) (Subscript[X, i]),  n ≠ 0}})

for consecutive points (n, (Subscript[S, n])),  n = 0, 1, 
      2, ...  
⟶a grah wich resembles the irregular motion of pollen grains.

1 - dimensional random walk

n = 20  0;
m = {0};
i = 1;
While[i < n, 
          m = Append[m, Extract[m, i] + 2*Random[Integer] - 1]; i++];
m1 = {0};
j = 1;
While[j < n, 
          m1 = Append[m1, 
              Extract[m1, j] + 2*Random[Integer] - 1]; j++];
m2 = {0};
k = 1;
While[k < n, 
          m2 = Append[m2, 
              Extract[m2, k] + 2*Random[Integer] - 1]; k++];
p1 = 
          ListPlot[m, PlotJoined -> True, DisplayFunction -> Identity, 
            PlotStyle -> Hue[ .43]];
p2 = 
          ListPlot[m1, PlotJoined -> True, DisplayFunction -> Identity, 
            PlotStyle -> Hue[ .95]];
p3 = 
          ListPlot[m2, PlotJoined -> True, DisplayFunction -> Identity, 
            PlotStyle -> Hue[ .75]];
Show[{p1, p2, p3},DisplayFunction -> $DisplayFunction,Axes -> True];

2 - dimensional random walk

n = 3000;                      
(Subscript[x, 0]) = 0;                           
(Subscript[y, 0]) = 0;                          
m = {0};
i = 1;
While[i < n + 1, 
          m = Append[m, Extract[m, i] + 2*Random[Integer] - 1]; i++];
j = 1;
s = {0};
While[j < n + 1, 
          s = Append[s, Extract[s, j] + 2*Random[Integer] - 1]; j++];
\(z = {}; \  (*\ 
        random\ walk\ with\ coordinate\ m\ \(an  d\)\ s\ *) 
j = 1;\)
While[j < n, 
          z = Append[z, {m[[j]], s[[j]]}]; j++];
ListPlot[z, PlotJoined -> True, PlotStyle -> Hue[ .44]]

* Since:

\(EY\_i = 0, \ i ∈ \(ℕ\^*\)\)

\(Var \((Y\_i)\) = \(\((E \((Y\_i - EY\_i)\))\)\^2 = 1\), \ 
      i ∈ \(ℕ\^*\)\)

obtain:

(Subscript[ES, n]) = 0

Var ((Subscript[S, n])) = (∑_(i = 1)^n)((Subscript[Y, i])) = n

and by the Central Limit Theorem (CLT)  we obtain that

((Subscript[S, [n]])/(Sqrt[n])) (Overscript[(n -> ∞), ⇀]) (Subscript[B, 1])

in distribution, where (Subscript[B, 1]) is random variable with a normal N(0, 1) distribution.

Construct extension for all times t ≥ 0

(Subscript[B, t]) = 
        (lim_(n -> ∞))((Subscript[S, [n t]])/(Sqrt[n ]))

1 - dimensional Brownian Motion

n = 10;
m = {0};
(Subscript[x, 0]) = 0;
t = 5;                
i = 1;
m = {0};
While[i < n*t + 1, 
          m = Append[m, Extract[m, i] + 2*Random[Integer] - 1]; i++];
 ListPlot[m,PlotJoined -> True,PlotStyle->Hue[.43],Axes -> True]; 
X[u_] := 
          (Subscript[x, 0]) + Extract[m, Max[IntegerPart[n*u], 1]]/Sqrt[n];
Plot[X[u],{u, 0, t},Compiled -> True,PlotStyle -> Hue[ .43],AspectRatio -> Automatic,Axes -> True];

2 - dimensional Brownian Motion

n = 300;                      
(Subscript[x, 0]) = 0;                           
(Subscript[y, 0]) = 0;                          
t = 5;       
m = {0};
i = 1;
While[i < n*t + 1,m=Append[m,Extract[m, i]+2*Random[Integer]-1];i++];
j = 1;
s = {0};
While[j < n*t + 1, 
          s = Append[s, Extract[s, j] + 2*Random[Integer] - 1]; j++];
z = {}; 
j = 1;
While[j < n*t, 
          z = Append[z, {m[[j]], s[[j]]}]; j++];
ListPlot[z, PlotJoined -> True, PlotStyle -> Hue[ .44]]
X[u_] := 
            (Subscript[x, 0]) + Extract[Extract[z, Max[IntegerPart[n*u], 1]], 1]/
                Sqrt[n];
  
Y[u_] := 
            (Subscript[y, 0]) + Extract[Extract[z, Max[IntegerPart[n*u], 1]], 2]/
                Sqrt[n];
  
ParametricPlot[{X[u], Y[u]}, {u, 0, t}, 
          PlotStyle -> Hue[ .44], AspectRatio -> Automatic];

Remark 
Random walk converges to Brownian motion (Donsker (1951))

Definition of Brownian Motion

Brownian motion is the unique process with the following properties:
(i) No memory ⟶ (Subscript[B, (Subscript[t, 1])])-(Subscript[B, (Subscript[t, 0])]),(Subscript[B, (Subscript[t, 2])]) - (Subscript[B, (Subscript[t, 1])]),(Subscript[B, (Subscript[t, 3])]) - (Subscript[B, (Subscript[t, 2])]),...  are independent
(ii) Invariance ⟶ the distribution of (Subscript[B, s + t])-(Subscript[B, t ])depends only on t
(iii) Continuity ⟶             t→(Subscript[B, t]) is continuous a . s and is nowhere differentiable a . s
(iv)  (Subscript[B, 0])=0, E ((Subscript[B, t])) = 0, Var ((Subscript[B, t])) = t

Multidimensional Brownian motion

(Subscript[B, t]^1),(Subscript[B, t]^2),(Subscript[B, t]^3),(Subscript[B, t]^4),… ⟶  independent 1-dimensional  Brownian motions
((Subscript[B, t]^1), (Subscript[B, t]^2), (Subscript[B, t]^3), ⋯,  
              (Subscript[B, t]^d))⟶  d-dimensional Brownian motion

Reflected Brownian motion RBM

* Reflected random walk converges to reflected Brownian motion 
(Stroock and Varadhan (1971) - (ℂ^2)  domains,Burdzy and Chen (2007) – uniform domains, not all domains)

Definition

(Subscript[X, t]) is a reflecting Brownian motion in D started at (Subscript[x, 0])∈D if it satisfies the stochastic differential equation

(Subscript[X, t]) = 
        (Subscript[X, 0]) + (Subscript[B, t]) + (1/2) (∫_0^t) (Subscript[λ, D]) ((Subscript[X, s])) (Subscript[dL, s]),

where

(Subscript[B, t])is a 2-dimensional  Brownian motion in D starting at0,
(Subscript[L, t])is a continuous nondecreasing process wich increases only when (Subscript[X, t])∈∂D
almost surly (Subscript[X, 0]) = 
        (Subscript[x, 0]) and (Subscript[X, t]) ∈ ∂D for all t ≥ 
            0

For 2 - dimensional

\(xincr = 1; \  (*\ the\ x\ increment\ *) 
\(yincr = 1;\)\)   
\(x\_0 = 0; \  (*\ 
        the\ initial\ x\  - \ coordinate\ of\ the\ random\ \(wal  k\)\ *) 

      y\_0 = 0; \  (*\ 
        the\ initial\ y\  - \ coordinate\ of\ the\ random\ walk\ *) 

      t\_0 = 3; \  (*\ 
        maximum\ value\ for\ the\ time\ variable\ t\ of\ the\ Brownian\ motion\ B\_t\ *) 
n = 300*t\_0;\)
i = 1;
m = {0};
While[i < n*t + 1, 
          m = Append[m, Extract[m, i] + 2*Random[Integer] - 1]; i++];
j = 1;
s = {0};
While[j < n*t + 1, 
          s = Append[s, Extract[s, j] + 2*Random[Integer] - 1]; j++];
z = {};
j = 1;
While[j < n*t, 
          z = Append[z, {m[[j]], s[[j]]}]; j++];
ListPlot[z, PlotJoined -> True, 
          PlotStyle -> Hue[ .44]]; 
localtime = 
          Rest[FoldList[Max, 0, 
              Table[-s[[i]] - (Subscript[y, 0])*Sqrt[n]*yincr, {i, n + 1}]]];
X[u_] := 
          (Subscript[x, 0]) + Extract[m, Max[IntegerPart[n*u], 1]]/Sqrt[n];
Y[u_] := 
          (Subscript[y, 0]) + Extract[s, Max[IntegerPart[n*u], 1]]/Sqrt[n];
Z[u_] := 
          (Subscript[y, 0]) + (Extract[s, Max[IntegerPart[n*u], 1]] + 
                  localtime[[1 + IntegerPart[n*u]]])/Sqrt[n];
\(ParametricPlot[{{X[u], Y[u]}, {X[u], Z[u]}}, {u, 0, 1}, 
          PlotStyle -> {{Hue[0.45], Thickness[0.006]}, {Hue[0.93], 
                Thickness[0.003]}}, 
          AspectRatio -> 
            Automatic, 
			PlotRange -> {{\(-1.2\), 1.2}, {\(-1.2\), 
                1.2}}];\)

The heat equation

Consider the heat equation in D = (ℝ^d)
Consider the problem of finding a (ℂ^(1, 
            2)) function

u = u (t, x) : [0, ∞)\[Cross]
          D → ℝ

satisfying the following initial value problem

{({{(1/2) (Subscript[Δ, x]) 
                        u (t, x) = ∂(u/∂)t (t, 
                          x),  (t, 
                        x) ∈ (0, +∞)\[Cross]
                        D, }, {u (0, 
                          x) = (Subscript[u, 0]) (x), x ∈ 
                      D,                                                   }})

where

(Subscript[Δ, x]) = ((∂^2)/∂)(Subscript[x, 1]^2) + … + ((∂^2)/∂)(Subscript[x, d]^2)

\(denotes\ the\ Laplacian\ operator\ in\ the\ variable\ x = \((x\_1, x\_2, …, x\_d)\) ∈ D\ and, \ 
      continuous\ function\ \(\(u\_0\) : 
          D\&_⟶ℝ\ represents\ the\ initial\ temperature\ distribution\ in\ \(\(D\)\(.\)\)\)\)

Remark.

u = u (t, x) : [0, ∞)\[Cross]
            D→ℝ          is the approximation solution of initial value problem, given by

u (t, x) = (E^x) ((Subscript[u, 0]) ((Subscript[B, t])))

with (Subscript[B, t])  -  d - 
          dimensional Brownian Motion starting at (Subscript[B, 0]) = x

* Problem for the heat equation in a bounded domain  D⊂(ℝ^d), with  Dirichlet conditions on the boundry (killed Brownian motion  - KBM) of D is

{({{(1/2) (Subscript[Δ, x]) 
                        u (t, x) = ∂(u/∂)t (t, 
                          x), (t, 
                        x) ∈ (0, +∞)\[Cross]
                        D,
u (t, x) = 
                      0, x ∈ ∂D,                                                           }, {u (0, 
                          x) = (Subscript[u, 0]) (x), x ∈ 
                      D,                                                    }})

Using probabilistic methods we can give explicit representation of solutions of heat equation with Dirichlet condition

u (t, x)=(E^x)((Subscript[u, 0])((Subscript[B, t])) (Subscript[1, t < (Subscript[τ, D])]))

where(Subscript[B, t])is a Brownian motion killed on the boundary of D starting at(Subscript[B, 0]) = x,
and(E^x)denotes the expectation of function (Subscript[u, 0])at(Subscript[B, t])with {t<(Subscript[τ, D])}

*Problem for the heat equation in a bounded domain D⊂(ℝ^d), with Neumann conditions on the boundary (reflected Brownian motion  - 
              RBM) of D is

{({{(1/2) (Subscript[Δ, x]) 
                        u (t, x) = ∂(u/∂)t (t, 
                          x), (t, 
                        x) ∈ (0, +∞)\[Cross]
                        D,
∂(u/∂)n (t, x) = 
                      0, x ∈ ∂D,                                                           }, {u (0, 
                          x) = (Subscript[u, 0]) (x), x ∈ 
                      D,                                                    }})

where \(\(u\_0\) : 
            D\&_⟶ℝ\ is\ a\ given\ continuous\ function\ \((the\ initial\ temperature\ distribution\ in\ D)\)\),
n = 
          n (x) denotes the inward normal unit vector field on ∂D.

Using probabilistic methods we can give explicit representation of solutions of heat equation with Neumann condition

u (t, x)=(E^x) (Subscript[u, 0])((Subscript[B, t]))

where(Subscript[B, t])is a reflected Brownian motion in D starting at(Subscript[B, 0]) = x,
and(E^x)denotes the expectation of function (Subscript[u, 0])at(Subscript[B, t])

*In (ℝ^2) we have

u (t, x)=(E^x) (Subscript[u, 0])((Subscript[B, t]))=(∫_((ℝ^2)))(Subscript[u, 0]) (y) 
                p (t, x, y) 
                dy = (∫_((ℝ^2)))(Subscript[u, 0]) (y) (1/(2  πt)) (E^(-((\[VerticalSeparator]x - y (\[VerticalSeparator] ^2))/(2  t)))) dy

p (t, x, 
            y) = (1/(2  πt)) (E^(-((\[VerticalSeparator]x - y (\[VerticalSeparator] ^2))/(2  
                      t)))) dy

with  x = ((Subscript[x, 1]), (Subscript[x, 2])),  y = ((Subscript[y, 1]), (Subscript[y, 2]))

* In H= upper half plane

u (t, x)=(E^x) (Subscript[u, 0])((Subscript[B, t]))=(∫_H) (Subscript[u, 0]) (y) (Subscript[p, H]) (t, x, y) 
              dy

where

* for Dirichlet Condition we have

\(\(p\_H\%D\) \((t, x, y)\) = 
        p \((t, x, y)\) - p \((t, x, y\&_)\), y\&_ = \((y\_1, \(-y\_2\))\)\)

* for Neumann Condition we have

\(\(p\_H\%N\) \((t, x, y)\) = 
        p \((t, x, y)\) + p \((t, x, y\&_)\), y\&_ = \((y\_1, \(-y\_2\))\)\)

Feynman - Kac formula:

u (t, x) = 
        E ((Subscript[u, 0]) ((Subscript[B, t]))) ≈ (((Subscript[u, 0]) ((Subscript[B, t]) ((Subscript[ω, 1]))) + (Subscript[u, 0]) ((Subscript[B, t]) ((Subscript[ω, 2]))) + … + (Subscript[u, 0]) ((Subscript[B, t]) ((Subscript[ω, n]))))/n)

Numerical results

For Neumann Condition

n = 1000;   
(Subscript[x, 0]) = 0;   
(Subscript[y, 0]) =  .3;   
t = 3;   
i = 1;
m = {0};
While[i < n*t + 1, 
        m = Append[m, 
            Extract[m, i] + 2*Random[Integer] - 1]; i++];
j = 
        1;
s = {0};
While[j < n*t + 1, 
        s = Append[s, 
            Extract[s, j] + 2*Random[Integer] - 1]; j++];
  
  

X[u_] := 
        (Subscript[x, 0]) + Extract[m, Max[IntegerPart[n*u], 1]]/Sqrt[n];
Y[
          u_] := 
        (Subscript[y, 0]) + Extract[s, Max[IntegerPart[n*u], 1]]/
            Sqrt[n];
\(ParametricPlot[{{X[u], Abs[Y[u]]}, {X[u], 
          Y[u]}}, {u, 0, 1}, PlotStyle -> {Hue[ .44], Hue[ .95]}, 
      AspectRatio -> Automatic]; \ \  (*\ RBM\ with\ green, \ 
      BM\ with\ red\ *) 
f[x_, y_] := 
      1\/\(x\^2 + y\^2\); \  (*\ u\_0 = f\ *) 
G = 0;\)
For[k = 1, k < n, G = G + f[X[t], Abs[Y[t]]], k++];
  
N[G/n]

n = 10;
(Subscript[x, 0]) = 0;
(Subscript[y, 0]) =  .3;
t = 5;
f[x_, y_] := (1/((x^2) + (y^2)));
X[x_, u_] := 
        Module[{i = 1, m = {0}}, 
          While[i < n*t + 1, 
            m = Append[m, Extract[m, i] + 2*Random[Integer] - 1]; i++]; 
          x + Extract[m, Max[IntegerPart[n*u], 1]]/
              Sqrt[n]];
Y[y_, u_] := 
      Module[{j = 1, s = {0}}, 
        While[j < n*t + 1, 
          s = Append[s, Extract[s, j] + 2*Random[Integer] - 1]; j++]; 
        y + Extract[s, Max[IntegerPart[n*u], 1]]/Sqrt[n]]

Density RBM  -  Neumann

\(\(p\_H\%N\) \((t, x, y)\) = 
        p \((t, x, y)\) + p \((t, x, y\&_)\), y\&_ = \((y\_1, \(-y\_2\))\)\)

u (t, x)=(E^x) (Subscript[u, 0])((Subscript[B, t]))=(∫_H) (Subscript[u, 0]) (y) (Subscript[p, H]) (t, x, y) 
              dy

p[t_, x1_, x2_, y1_, y2_] = 
        (1/(2  π t))*
            Exp[-(((((x1 - y1))^2) + (((x2 - y2))^2))/(2  t))] + 
          (1/(2  π t))*
            Exp[-(((((x1 - y1))^2) + (((x2 + y2))^2))/(2  
                      t))];
u[x1_, x2_, t_] = 
      Integrate[
        f[y1, y2] 
          p[t, x1, x2, y1, y2], {y1, -∞, ∞}, {y2, 
          0, ∞}]

v[x1_, x2_, t_] := 
        Module[{G, k}, G = 0; 
          For[k = 1, k < n, G = G + f[X[x1, t], Abs[Y[x2, t]]], k++]; 
          N[G/n]];

N[v[0, 0.3, 5]]

Remark
in the same way for heat equation with  Dirichlet conditions

References
[1] Pascu, M., "Brownian motion and application",  Transilvania University, Brasov, 2006
[2] Oksendal, B., "Stochastic differential equations: An introduction applications", six edition, Springer, 2003