Gallery of proved programs

dirichlet: discretization of the 1D acoustic wave equation

This example was developed among the FOST project. It is a finite difference numerical scheme for the resolution of the one-dimensional acoustic wave equation. A reasonable bound on the rounding error requires a complex property that expresses each rounding error. The bound on the method error of this scheme has also been formally certified in Coq. For this proof, we also require the Coquelicot Coq library.


/*@ axiomatic dirichlet_maths {
  @
  @ logic real c;
  @ logic real p0(real x);
  @ logic real psol(real x, real t);

  @ axiom c_pos: 0 < c;

  @ logic real psol_1(real x, real t);
  @ axiom psol_1_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dx;
  @   0 < eps ==> \abs(dx) < C ==>
  @   \abs((psol(x + dx, t) - psol(x, t)) / dx - psol_1(x, t)) < eps;

  @ logic real psol_11(real x, real t);
  @ axiom psol_11_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dx;
  @   0 < eps ==> \abs(dx) < C ==>
  @   \abs((psol_1(x + dx, t) - psol_1(x, t)) / dx - psol_11(x, t)) < eps;

  @ logic real psol_2(real x, real t);
  @ axiom psol_2_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dt;
  @   0 < eps ==> \abs(dt) < C ==>
  @   \abs((psol(x, t + dt) - psol(x, t)) / dt - psol_2(x, t)) < eps;

  @ logic real psol_22(real x, real t);
  @ axiom psol_22_def:
  @   \forall real x; \forall real t;
  @   \forall real eps; \exists real C; 0 < C && \forall real dt;
  @   0 < eps ==> \abs(dt) < C ==>
  @   \abs((psol_2(x, t + dt) - psol_2(x, t)) / dt - psol_22(x, t)) < eps;

  @ axiom wave_eq_0: \forall real x; 0 <= x <= 1 ==> psol(x, 0) == p0(x);
  @ axiom wave_eq_1: \forall real x; 0 <= x <= 1 ==> psol_2(x, 0) == 0;
  @ axiom wave_eq_2:
  @   \forall real x; \forall real t;
  @   0 <= x <= 1 ==> psol_22(x, t) - c * c * psol_11(x, t) == 0;
  @ axiom wave_eq_dirichlet_1: \forall real t; psol(0, t) == 0;
  @ axiom wave_eq_dirichlet_2: \forall real t; psol(1, t) == 0;

  @ logic real psol_Taylor_3(real x, real t, real dx, real dt);
  @ logic real psol_Taylor_4(real x, real t, real dx, real dt);

  @ logic real alpha_3; logic real C_3;
  @ logic real alpha_4; logic real C_4;

  @ axiom psol_suff_regular_3:
  @   0 < alpha_3 && 0 < C_3 &&
  @   \forall real x; \forall real t; \forall real dx; \forall real dt;
  @   0 <= x <= 1 ==> \sqrt(dx * dx + dt * dt) <= alpha_3 ==>
  @   \abs(psol(x + dx, t + dt) - psol_Taylor_3(x, t, dx, dt)) <=
  @     C_3 * \abs(\pow(\sqrt(dx * dx + dt * dt), 3));

  @ axiom psol_suff_regular_4:
  @   0 < alpha_4 && 0 < C_4 &&
  @   \forall real x; \forall real t; \forall real dx; \forall real dt;
  @   0 <= x <= 1 ==> \sqrt(dx * dx + dt * dt) <= alpha_4 ==>
  @   \abs(psol(x + dx, t + dt) - psol_Taylor_4(x, t, dx, dt)) <=
  @     C_4 * \abs(\pow(\sqrt(dx * dx + dt * dt), 4));

  @ axiom psol_le:
  @   \forall real x; \forall real t;
  @   0 <= x <= 1 ==> 0 <= t ==> \abs(psol(x, t)) <= 1;

  @ logic real T_max;
  @ axiom T_max_pos: 0 < T_max;

  @ logic real C_conv; logic real alpha_conv;
  @ lemma alpha_conv_pos: 0 < alpha_conv;
  @
  @ } */

/*@ axiomatic dirichlet_prog {
  @
  @  predicate analytic_error{L}
  @    (double **p, integer ni, integer i, integer k, double a, double dt)
  @    reads p[..][..];
  @
  @  lemma analytic_error_le{L}:
  @    \forall double **p; \forall integer ni; \forall integer i;
  @    \forall integer nk; \forall integer k;
  @    \forall double a; \forall double dt;
  @    0 < ni ==> 0 <= i <= ni ==> 0 <= k ==>
  @    0 < \exact(dt) ==>
  @    analytic_error(p, ni, i, k, a, dt) ==>
  @    \sqrt(1. / (ni * ni) + \exact(dt) * \exact(dt)) < alpha_conv ==>
  @    k <= nk ==> nk <= 7598581 ==> nk * \exact(dt) <= T_max ==>
  @    \exact(dt) * ni * c <= 1 - 0x1.p-50 ==>
  @    \forall integer i1; \forall integer k1;
  @    0 <= i1 <= ni ==> 0 <= k1 < k ==>
  @    \abs(p[i1][k1]) <= 2;
  @
  @  predicate separated_matrix{L}(double **p, integer leni) =
  @    \forall integer i; \forall integer j;
  @    0 <= i < leni ==> 0 <= j < leni ==> i != j ==>
  @    \base_addr(p[i]) != \base_addr(p[j]);
  @
  @ logic real sqr_norm_dx_conv_err{L}
  @   (double **p, real dx, real dt, integer ni, integer i, integer k)
  @   reads p[..][..];
  @ logic real sqr(real x) = x * x;
  @ lemma sqr_norm_dx_conv_err_0{L}:
  @   \forall double **p; \forall real dx; \forall real dt;
  @   \forall integer ni; \forall integer k;
  @   sqr_norm_dx_conv_err(p, dx, dt, ni, 0, k) == 0;
  @ lemma sqr_norm_dx_conv_err_succ{L}:
  @   \forall double **p; \forall real dx; \forall real dt;
  @   \forall integer ni; \forall integer i; \forall integer k;
  @   0 <= i ==>
  @   sqr_norm_dx_conv_err(p, dx, dt, ni, i + 1, k) ==
  @     sqr_norm_dx_conv_err(p, dx, dt, ni, i, k) +
  @     dx * sqr(psol(1. * i / ni, k * dt) - \exact(p[i][k]));
  @ logic real norm_dx_conv_err{L}
  @   (double **p, real dt, integer ni, integer k) =
  @   \sqrt(sqr_norm_dx_conv_err(p, 1. / ni, dt, ni, ni, k));
  @
  @ } */

/*@ requires leni >= 1 && lenj >= 1;
  @ ensures
  @   \valid_range(\result, 0, leni - 1) &&
  @   (\forall integer i; 0 <= i < leni ==>
  @    \valid_range(\result[i], 0, lenj - 1)) &&
  @   separated_matrix(\result, leni);
  @ */
double **array2d_alloc(int leni, int lenj);

/*@ requires (l != 0) && \round_error(x) <= 5./2*0x1.p-53;
  @ ensures
  @   \round_error(\result) <= 14 * 0x1.p-52 &&
  @   \exact(\result) == p0(\exact(x));
  @ */
double p_zero(double xs, double l, double x);

/*@ requires
  @   ni >= 2 && nk >= 2 && l != 0 &&
  @   dt > 0. && \exact(dt) > 0. &&
  @   \exact(v) == c && \exact(v) == v &&
  @   0x1.p-1000 <= \exact(dt) &&
  @   ni <= 2147483646 && nk <= 7598581 &&
  @   nk * \exact(dt) <= T_max &&
  @   \abs(\exact(dt) - dt) / dt <= 0x1.p-51 &&
  @   0x1.p-500 <= \exact(dt) * ni * c <= 1 - 0x1.p-50 &&
  @   \sqrt(1. / (ni * ni) + \exact(dt) * \exact(dt)) < alpha_conv;
  @
  @ ensures
  @   \forall integer i; \forall integer k;
  @   0 <= i <= ni ==> 0 <= k <= nk ==>
  @   \round_error(\result[i][k]) <= 78. / 2 * 0x1.p-52 * (k + 1) * (k + 2);
  @
  @ ensures
  @   \forall integer k; 0 <= k <= nk ==>
  @   norm_dx_conv_err(\result, \exact(dt), ni, k) <=
  @     C_conv * (1. / (ni * ni) + \exact(dt) * \exact(dt));
  @ */
double **forward_prop(int ni, int nk, double dt, double v, double xs, double l) {

  /* Output variable. */
  double **p;

  /* Local variables. */
  int i, k;
  double a1, a, dp, dx;

  dx = 1./ni;
  /*@ assert
    @   dx > 0. && dx <= 0.5 &&
    @   \abs(\exact(dx) - dx) / dx <= 0x1.p-53;
    @ */

  /* Compute the constant coefficient of the stiffness matrix. */
  a1 = dt/dx*v;
  a = a1*a1;
  /*@ assert
    @   0 <= a <= 1 &&
    @   0 < \exact(a) <= 1 &&
    @   \round_error(a) <= 0x1.p-49;
    @ */

  /* Allocate space-time variable for the discrete solution. */
  p = array2d_alloc(ni+1, nk+1);

  /* First initial condition and boundary conditions. */
  /* Left boundary. */
  p[0][0] = 0.;
  /* Time iteration -1 = space loop. */
  /*@ loop invariant
    @   1 <= i <= ni &&
    @   analytic_error(p, ni, i - 1, 0, a, dt);
    @ loop variant ni - i; */
  for (i=1; i<ni; i++) {
    p[i][0] = p_zero(xs, l, i*dx);
  }
  /* Right boundary. */
  p[ni][0] = 0.;
  /*@ assert analytic_error(p, ni, ni, 0, a, dt); */

  /* Second initial condition (with p_one=0) and boundary conditions. */
  /* Left boundary. */
  p[0][1] = 0.;
  /* Time iteration 0 = space loop. */
  /*@ loop invariant
    @   1 <= i <= ni &&
    @   analytic_error(p, ni, i - 1, 1, a, dt);
    @ loop variant ni - i; */
  for (i=1; i<ni; i++) {
    /*@ assert \abs(p[i-1][0]) <= 2; */
    /*@ assert \abs(p[i][0]) <= 2; */
    /*@ assert \abs(p[i+1][0]) <= 2; */
    dp = p[i+1][0] - 2.*p[i][0] + p[i-1][0];
    p[i][1] = p[i][0] + 0.5*a*dp;
  }
  /* Right boundary. */
  p[ni][1] = 0.;
  /*@ assert analytic_error(p, ni, ni, 1, a, dt); */

  /* Evolution problem and boundary conditions. */
  /* Propagation = time loop. */
  /*@ loop invariant
    @   1 <= k <= nk &&
    @   analytic_error(p, ni, ni, k, a, dt);
    @ loop variant nk - k; */
  for (k=1; k<nk; k++) {
    /* Left boundary. */
    p[0][k+1] = 0.;
    /* Time iteration k = space loop. */
    /*@ loop invariant
      @   1 <= i <= ni &&
      @   analytic_error(p, ni, i - 1, k + 1, a, dt);
      @ loop variant ni - i; */
    for (i=1; i<ni; i++) {
      /*@ assert \abs(p[i-1][k]) <= 2; */
      /*@ assert \abs(p[i][k]) <= 2; */
      /*@ assert \abs(p[i+1][k]) <= 2; */
      /*@ assert \abs(p[i][k-1]) <= 2; */
      dp = p[i+1][k] - 2.*p[i][k] + p[i-1][k];
      p[i][k+1] = 2.*p[i][k] - p[i][k-1] + a*dp;
    }
    /* Right boundary. */
    p[ni][k+1] = 0.;
    /*@ assert analytic_error(p, ni, ni, k + 1, a, dt); */
  }

  return p;

}