FreeFEM Documentation on GitHub

stars - forks

Evolution problems

FreeFEM also solves evolution problems such as the heat equation:

(60)\[\begin{split}\begin{array}{rcll} \frac{\p u}{\p t}-\mu\Delta u &=& f & \textrm{ in }\Omega\times ]0,T[\\ u(\mathbf{x},0) &=& u_0(\mathbf{x}) & \textrm{ in }\Omega\\ \left(\p u/\p n\right)(\mathbf{x},t) &=& 0 & \textrm{ on }\p\Omega\times ]0,T[ \end{array}\end{split}\]

with a positive viscosity coefficient \(\mu\) and homogeneous Neumann boundary conditions.

We solve (60) by FEM in space and finite differences in time.

We use the definition of the partial derivative of the solution in the time derivative:

\[\frac{\p u}{\p t}(x,y,t) = \lim_{\tau \to 0}\frac{u(x,y,t)-u(x,y,t-\tau )}{\tau }\]

which indicates that \(u^m(x,y)=u(x,y,m\tau )\) will satisfy approximatively:

\[\frac{\p u}{\p t}(x,y,m\tau )\simeq \frac{u^m(x,y)-u^{m-1}(x,y)}{\tau }\]

The time discretization of heat equation (60) is as follows, \(\forall m=0,\cdots,[T/\tau ]\):

\[\begin{split}\begin{array}{rcll} \frac{u^{m+1}-u^{m}}{\tau }-\mu\Delta u^{m+1} &=& f^{m+1} & \textrm{ in }\Omega\\ u^0(\mathbf{x}) &=& u_0(\mathbf{x}) & \textrm{ in }\Omega\\ \p u^{m+1}/\p n(\mathbf{x}) &=& 0 & \textrm{ on }\p\Omega \end{array}\end{split}\]

which is so-called backward Euler method for (60).

To obtain the variational formulation, multiply with the test function \(v\) both sides of the equation:

\[\int_{\Omega}\{u^{m+1}v-\tau \Delta u^{m+1}v\}=\int_{\Omega}\{u^m+\tau f^{m+1}\}v\]

By the divergence theorem, we have:

\[\int_{\Omega}\{u^{m+1}v+\tau\nabla u^{m+1}\cdot \nabla v\} -\int_{\p\Omega} \tau \left( \p u^{m+1}/\p n\right) v =\int_{\Omega }\{u^mv+\tau f^{m+1}v\}\]

By the boundary condition \(\p u^{m+1}/\p n=0\), it follows that:

(61)\[\int_{\Omega} \{u^{m+1}v+\tau \nabla u^{m+1}\cdot \nabla v\} -\int_{\Omega }\{u^mv+\tau f^{m+1}v\} =0\]

Using the identity just above, we can calculate the finite element approximation \(u_h^m\) of \(u^m\) in a step-by-step manner with respect to \(t\).

Tip

Example

We now solve the following example with the exact solution \(u(x,y,t)=tx^4\), \(\Omega = ]0,1[^2\).

\[\begin{split}\begin{array}{rcll} \frac{{\p u}}{{\p t}} - \mu \Delta u &=& x^4 - \mu 12tx^2 & \textrm{ in }\Omega\times ]0,3[\\ u(x,y,0) &=& 0 & \textrm{ on }\Omega\\ \left. u \right|_{\p\Omega} &=& t*x^4 \end{array}\end{split}\]
 1// Parameters
 2real dt = 0.1;
 3real mu = 0.01;
 4
 5// Mesh
 6mesh Th = square(16, 16);
 7
 8// Fespace
 9fespace Vh(Th, P1);
10Vh u, v, uu, f, g;
11
12// Problem
13problem dHeat (u, v)
14    = int2d(Th)(
15          u*v
16        + dt*mu*(dx(u)*dx(v) + dy(u)*dy(v))
17    )
18    + int2d(Th)(
19        - uu*v
20        - dt*f*v
21    )
22    + on(1, 2, 3, 4, u=g)
23    ;
24
25// Time loop
26real t = 0;
27uu = 0;
28for (int m = 0; m <= 3/dt; m++){
29    // Update
30    t = t+dt;
31    f = x^4 - mu*t*12*x^2;
32    g = t*x^4;
33    uu = u;
34
35    // Solve
36    dHeat;
37
38    // Plot
39    plot(u, wait=true);
40    cout << "t=" << t << " - L^2-Error=" << sqrt(int2d(Th)((u-t*x^4)^2)) << endl;
41}

In the last statement, the \(L^2\)-error \(\left(\int_{\Omega}\left| u-tx^4\right|^2\right)^{1/2}\) is calculated at \(t=m\tau, \tau =0.1\). At \(t=0.1\), the error is 0.000213269. The errors increase with \(m\) and 0.00628589 at \(t=3\).

The iteration of the backward Euler (61) is made by for loop.

Note

The stiffness matrix in the loop is used over and over again. FreeFEM support reuses of stiffness matrix.

Mathematical Theory on Time Difference Approximations.

In this section, we show the advantage of implicit schemes. Let \(V, H\) be separable Hilbert space and \(V\) is dense in \(H\). Let \(a\) be a continuous bilinear form over \(V \times V\) with coercivity and symmetry.

Then \(\sqrt{a(v,v)}\) become equivalent to the norm \(\| v\|\) of \(V\).

Problem Ev(f,Omega): For a given \(f\in L^2(0,T;V'),\, u^0\in H\)

\[\begin{split}\begin{array}{rcl} \frac{d}{dt}(u(t),v)+a(u(t),v)&=&( f(t),v)\qquad \forall v\in V,\quad a.e. \, t\in [0,T]\\ u(0)&=&u^0\nonumber \end{array}\end{split}\]

where \(V'\) is the dual space of \(V\).

Then, there is an unique solution \(u\in L^{\infty}(0,T;H)\cap L^2(0,T;V)\).

Let us denote the time step by \(\tau>0\), \(N_T=[T/\tau]\). For the discretization, we put \(u^n = u(n\tau)\) and consider the time difference for each \(\theta\in [0,1]\)

(62)\[\begin{split}\begin{array}{rcl} \frac{1}{\tau}\left( u_h^{n+1}-u_h^n,\phi_i\right) +a\left( u_h^{n+\theta},\phi_i\right)&=&\langle f^{n+\theta},\phi_i\rangle\\ i=1,&\cdots&, m,\quad n=0,\cdots, N_T\nonumber\\ u_h^{n+\theta}&=&\theta u_h^{n+1}+(1-\theta)u_h^n,\\ f^{n+\theta}&=&\theta f^{n+1}+(1-\theta)f^n\nonumber \end{array}\end{split}\]

Formula (62) is the forward Euler scheme if \(\theta=0\), Crank-Nicolson scheme if \(\theta=1/2\), the backward Euler scheme if \(\theta=1\).

Unknown vectors \(u^n=(u_h^1,\cdots,u_h^M)^T\) in

\[u_h^n(x)=u^n_1\phi_1(x)+\cdots+u^n_m\phi_m(x),\quad u^n_1,\cdots,u^n_m\in \R\]

are obtained from solving the matrix

(63)\[\begin{split}(M+\theta\tau A)u^{n+1}=\{M-(1-\theta)\tau A\}u^n +\tau\left\{\theta f^{n+1}+(1-\theta)f^n\right\}\\ M=(m_{ij}),\quad m_{ij}=(\phi_j,\phi_i),\qquad A=(a_{ij}),\quad a_{ij}=a(\phi_j,\phi_i)\nonumber\end{split}\]

Refer [TABATA1994], pp.70–75 for solvability of (63). The stability of (63) is in [TABATA1994], Theorem 2.13:

Let \(\{\mathcal{T}_h\}_{h\downarrow 0}\) be regular triangulations (see Regular Triangulation). Then there is a number \(c_0>0\) independent of \(h\) such that,

\[\begin{split}|u_h^n|^2\le \left\{ \begin{array}{lr} \frac{1}{\delta}\left\{ |u^0_h|^2+\tau \sum_{k=0}^{n-1}\|f^{k+\theta}\|^2_{V_h'} \right\}&\theta\in [0,1/2)\\ |u^0_h|^2+\tau \sum_{k=0}^{n-1}\|f^{k+\theta}\|^2_{V_h'}&\theta\in [1/2,1] \end{array} \right.\end{split}\]

if the following are satisfied:

  1. When \(\theta\in [0,1/2)\), then we can take a time step \(\tau\) in such a way that

    \[\tau <\frac{2(1-\delta)}{(1-2\theta)c_0^2}h^2\]

    for arbitrary \(\delta\in (0,1)\).

  2. When \(1/2\leq \theta\leq 1\), we can take \(\tau\) arbitrary.

Tip

Example

 1// Parameters
 2real tau = 0.1; real
 3theta = 0.;
 4
 5// Mesh
 6mesh Th = square(12, 12);
 7
 8// Fespace
 9fespace Vh(Th, P1);
10Vh u, v, oldU;
11Vh f1, f0;
12
13fespace Ph(Th, P0);
14Ph h = hTriangle; // mesh sizes for each triangle
15
16// Function
17func real f (real t){
18    return x^2*(x-1)^2 + t*(-2 + 12*x - 11*x^2 - 2*x^3 + x^4);
19}
20
21// File
22ofstream out("err02.csv"); //file to store calculations
23out << "mesh size = " << h[].max << ", time step = " << tau << endl;
24for (int n = 0; n < 5/tau; n++)
25    out << n*tau << ",";
26out << endl;
27
28// Problem
29problem aTau (u, v)
30    = int2d(Th)(
31          u*v
32        + theta*tau*(dx(u)*dx(v) + dy(u)*dy(v) + u*v)
33    )
34    - int2d(Th)(
35          oldU*v
36        - (1-theta)*tau*(dx(oldU)*dx(v) + dy(oldU)*dy(v) + oldU*v)
37    )
38    - int2d(Th)(
39          tau*(theta*f1 + (1-theta)*f0)*v
40    )
41    ;
42
43// Theta loop
44while (theta <= 1.0){
45    real t = 0;
46    real T = 3;
47    oldU = 0;
48    out << theta << ",";
49    for (int n = 0; n < T/tau; n++){
50        // Update
51        t = t + tau;
52        f0 = f(n*tau);
53        f1 = f((n+1)*tau);
54
55        // Solve
56        aTau;
57        oldU = u;
58
59        // Plot
60        plot(u);
61
62        // Error
63        Vh uex = t*x^2*(1-x)^2; //exact solution = tx^2(1-x)^2
64        Vh err = u - uex; // err = FE-sol - exact
65        out << abs(err[].max)/abs(uex[].max) << ",";
66    }
67    out << endl;
68    theta = theta + 0.1;
69}
../_images/EvolutionProblems_TimeDifference.png

Fig. 161 \(\max_{x\in\Omega}\vert u_h^n(\theta)-u_{ex}(n\tau)\vert\max_{x\in\Omega}\vert u_{ex}(n\tau)\vert at n=0,1,\cdots,29\)

We can see in Fig. 161 that \(u_h^n(\theta)\) become unstable at \(\theta=0.4\), and figures are omitted in the case \(\theta<0.4\).

Convection

The hyperbolic equation

(64)\[\p_t u +\mathbf{\alpha} \cdot \nabla u=f;\ \textrm{ for a vector-valued function }\mathbf{\alpha}\]

appears frequently in scientific problems, for example in the Navier-Stokes equations, in the Convection-Diffusion equation, etc.

In the case of 1-dimensional space, we can easily find the general solution \((x,t)\mapsto u(x,t)=u^0(x-\alpha t)\) of the following equation, if \(\alpha\) is constant,

(65)\[\begin{split}\begin{array}{rcl} \p_t u +\alpha\p_x u &=& 0\\ u(x,0) &=& u^0(x), \end{array}\end{split}\]

because \(\p_t u +\alpha\p_x u=-\alpha\dot{u}^0+a\dot{u}^0=0\), where \(\dot{u}^0=du^0(x)/dx\).

Even if \(\alpha\) is not constant, the construction works on similar principles. One begins with the ordinary differential equation (with the convention that \(\alpha\) is prolonged by zero apart from \((0,L)\times (0,T)\)):

\[\dot{X}(\tau )=+\alpha(X(\tau ),\tau ),\ \tau \in (0,t)\quad X(t)=x\]

In this equation \(\tau\) is the variable and \(x,t\) are parameters, and we denote the solution by \(X_{x,t}(\tau )\). Then it is noticed that \((x,t)\rightarrow v(X(\tau),\tau)\) in \(\tau=t\) satisfies the equation

\[\p _{t}v+\alpha\p _{x}v=\p _{t}X\dot{v}+a\p _{x}X\dot{v}% =0\]

and by the definition \(\p _{t}X=\dot{X}=+\alpha\) and \(\p_{x}X=\p _{x}x\) in \(\tau=t\), because if \(\tau =t\) we have \(X(\tau )=x\).

The general solution of (65) is thus the value of the boundary condition in \(X_{x, t}(0)\), that is to say \(u(x,t)=u^{0}(X_{x,t}(0))\) where \(X_{x,t}(0)\) is on the \(x\) axis, \(u(x,t)=u^{0}(X_{x,t}(0))\) if \(X_{x,t}(0)\) is on the axis of \(t\).

In higher dimension \(\Omega \subset R^{d},~d=2,3\), the equation for the convection is written

\[\p _{t}u+\mathbf{\alpha}\cdot \nabla u=0\hbox{ in }\Omega \times (0,T)\]

where \(\mathbf{a}(x,t)\in \R^{d}\).

FreeFEM implements the Characteristic-Galerkin method for convection operators. Recall that the equation (64) can be discretized as

\[\frac{Du}{Dt} = f\;\;\textrm{i.e. }\frac{du}{dt}\left( {X(t),t} \right) = f\left(X( t ),t \right)\textrm{ where }\frac{dX}{dt}( t ) = \mathbf{\alpha}( {X(t),t})\]

where \(D\) is the total derivative operator. So a good scheme is one step of backward convection by the method of Characteristics-Galerkin

(66)\[\frac{1}{{\tau }}\left(u^{m + 1}(x) - u^m(X^m(x))\right) = f^{m+1} (x)\]

where \(X^m (x)\) is an approximation of the solution at \(t = m\tau\) of the ordinary differential equation

\[\frac{d\mathbf{X}}{dt}(t) = \mathbf{\alpha}^m(\mathbf{X}(t)), \mathbf{X}((m + 1)\tau) = x.\]

where \(\mathbf{\alpha}^m(x)=(\alpha_1(x,m\tau ),\alpha_2(x,m\tau))\). Because, by Taylor’s expansion, we have

(67)\[\begin{split}\begin{array}{rcl} u^m(\mathbf{X}(m\tau ))&=& u^m(\mathbf{X}((m+1)\tau )) - \tau \sum_{i=1}^d \frac{\p u^m}{\p x_i}(\mathbf{X}((m+1)\tau )) \frac{\p X_i}{\p t}((m+1)\tau ) +o(\tau )\nonumber\\ &=&u^m(x)-\tau \mathbf{\alpha}^m(x)\cdot \nabla u^m(x)+o(\tau ) \end{array}\end{split}\]

where \(X_i(t)\) are the i-th component of \(\mathbf{X}(t)\), \(u^m(x)=u(x,m\tau )\) and we used the chain rule and \(x=\mathbf{X}((m+1)\tau )\). From (67), it follows that

\[u^m(X^m(x))=u^m(x)-\tau \mathbf{\alpha}^m(x)\cdot \nabla u^m(x)+o(\tau )\]

Also we apply Taylor’s expansion for \(t \rightarrow u^m(x-\mathbf{\alpha}^m(x)t),0\le t\le \tau\), then

\[u^m(x-\mathbf{\alpha}\tau )=u^m(x)-\tau \mathbf{\alpha}^m(x)\cdot \nabla u^m(x)+o(\tau ).\]

Putting

convect\(\left( {\mathbf{\alpha},-\tau ,u^m } \right)\approx u^m \left(x - \mathbf{\alpha}^m\tau \right)\)

we can get the approximation

\(u^m \left( {X^m( x )} \right) \approx\) convect \(\left( {[a_1^m ,a_2^m],-\tau ,u^m } \right)\) by \(X^m \approx x \mapsto x- \tau [a_1^m(x) ,a_2^m(x)]\)

A classical convection problem is that of the “rotating bell” (quoted from [LUCQUIN1998], p.16).

Let \(\Omega\) be the unit disk centered at 0, with its center rotating with speed \(\alpha_1 = y,\, \alpha_2 = -x\). We consider the problem (64) with \(f=0\) and the initial condition \(u(x,0)=u^0(x)\), that is, from (66)

\(u^{m + 1}(x) = u^m(X^m(x))\approx\) convect\((\mathbf{\alpha},-\tau ,u^m)\)

The exact solution is \(u(x, t) = u(\mathbf{X}(t))\) where \(\mathbf{X}\) equals \(x\) rotated around the origin by an angle \(\theta = -t\) (rotate in clockwise). So, if \(u^0\) in a 3D perspective looks like a bell, then \(u\) will have exactly the same shape, but rotated by the same amount. The program consists in solving the equation until \(T = 2\pi\), that is for a full revolution and to compare the final solution with the initial one; they should be equal.

Tip

Convect

 1// Parameters
 2real dt = 0.17;
 3
 4// Mesh
 5border C(t=0, 2*pi){x=cos(t); y=sin(t);}
 6mesh Th = buildmesh(C(70));
 7
 8// Fespace
 9fespace Vh(Th, P1);
10Vh u0;
11Vh a1 = -y, a2 = x; //rotation velocity
12Vh u;
13
14// Initialization
15u = exp(-10*((x-0.3)^2 +(y-0.3)^2));
16
17// Time loop
18real t = 0.;
19for (int m = 0; m < 2*pi/dt; m++){
20    // Update
21    t += dt;
22    u0 = u;
23
24    // Convect
25    u = convect([a1, a2], -dt, u0); //u^{m+1}=u^m(X^m(x))
26
27    // Plot
28    plot(u, cmm=" t="+t+", min="+u[].min+", max="+u[].max);
29}

Note

The scheme convect is unconditionally stable, then the bell become lower and lower (the maximum of \(u^{37}\) is \(0.406\) as shown in Fig. 162.

EvolutionProblem_Convect

Fig. 162 \(u^0=e^{-10((x-0.3)^2 +(y-0.3)^2)}\)

EvolutionProblem_Convect2

Fig. 163 The bell at \(t=6.29\)

2D Black-Scholes equation for an European Put option

In mathematical finance, an option on two assets is modeled by a Black-Scholes equations in two space variables, (see for example [WILMOTT1995] or [ACHDOU2005]).

\[\begin{split}\begin{array}{rcl} \p _t u &+& \frac{{\left( {\sigma _1 x } \right)^2 }}{2}\frac{{\p ^2 u}}{{\p x^2 }} + \frac{{\left( {\sigma _2 y } \right)^2 }}{2}\frac{{\p ^2 u}}{{\p y^2 }} \\ &&{\rm{ }} + \rho x y \frac{{\p ^2 u}}{{\p x \p y }} + rS_1 \frac{{\p u}}{{\p x }} + rS_2 \frac{{\p u}}{{\p y }} - rP = 0 \nonumber \end{array}\end{split}\]

which is to be integrated in \(\left( {0,T} \right) \times \R^ + \times \R^ +\) subject to, in the case of a put

\[u\left( {x , y ,T} \right) = \left( {K - \max \left( {x ,y } \right)} \right)^+\]

Boundary conditions for this problem may not be so easy to device. As in the one dimensional case the PDE contains boundary conditions on the axis \(x_1 = 0\) and on the axis \(x_2 = 0\), namely two one dimensional Black-Scholes equations driven respectively by the data \(u\left( {0, + \infty ,T} \right)\) and \(u\left( { + \infty ,0,T} \right)\). These will be automatically accounted for because they are embedded in the PDE. So if we do nothing in the variational form (i.e. if we take a Neumann boundary condition at these two axis in the strong form) there will be no disturbance to these. At infinity in one of the variable, as in 1D, it makes sense to impose \(u=0\). We take

\[\sigma _1 = 0.3,\;\;\sigma _2 = 0.3,\;\;\rho = 0.3,\;\;r = 0.05,\;\;K = 40,\;\;T = 0.5\]

An implicit Euler scheme is used and a mesh adaptation is done every 10 time steps. To have an unconditionally stable scheme, the first order terms are treated by the Characteristic Galerkin method, which, roughly, approximates

\[\frac{{\p u}}{{\p t}} + a_1 \frac{{\p u}}{{\p x}} + a_2 \frac{{\p u}}{{\p y}} \approx \frac{1}{{\tau }}\left( {u^{n + 1} \left( x \right) - u^n \left( {x - \mathbf{\alpha}\tau } \right)} \right)\]

Tip

Black-Scholes

 1// Parameters
 2int m = 30; int L = 80; int LL = 80; int j = 100; real sigx = 0.3; real sigy = 0.3; real rho = 0.3; real r = 0.05; real K = 40; real dt = 0.01;
 3
 4// Mesh
 5mesh th = square(m, m, [L*x, LL*y]);
 6
 7// Fespace
 8fespace Vh(th, P1);
 9Vh u = max(K-max(x,y),0.);
10Vh xveloc, yveloc, v, uold;
11
12// Time loop
13for (int n = 0; n*dt <= 1.0; n++){
14    // Mesh adaptation
15    if (j > 20){
16        th = adaptmesh(th, u, verbosity=1, abserror=1, nbjacoby=2,
17        err=0.001, nbvx=5000, omega=1.8, ratio=1.8, nbsmooth=3,
18        splitpbedge=1, maxsubdiv=5, rescaling=1);
19        j = 0;
20        xveloc = -x*r + x*sigx^2 + x*rho*sigx*sigy/2;
21        yveloc = -y*r + y*sigy^2 + y*rho*sigx*sigy/2;
22        u = u;
23    }
24
25    // Update
26    uold = u;
27
28    // Solve
29    solve eq1(u, v, init=j, solver=LU)
30        = int2d(th)(
31              u*v*(r+1/dt)
32            + dx(u)*dx(v)*(x*sigx)^2/2
33            + dy(u)*dy(v)*(y*sigy)^2/2
34            + (dy(u)*dx(v) + dx(u)*dy(v))*rho*sigx*sigy*x*y/2
35        )
36        - int2d(th)(
37              v*convect([xveloc, yveloc], dt, uold)/dt
38        )
39        + on(2, 3, u=0)
40        ;
41
42    // Update
43    j = j+1;
44};
45
46// Plot
47plot(u, wait=true, value=true);

Results are shown on Fig. 164 and Fig. 165.

EvolutionProblems_BlackSholes

Fig. 164 The adapted triangulation

EvolutionProblems_BlackSholes2

Fig. 165 The level line of the European basquet put option

Table of content