Sylvie Boldo

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.

#pragma JessieIntegerModel(math)

/*@ axiomatic dirichlet {
  @  predicate separated_matrix(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]);
  @
  @  predicate analytic_error{L}
  @          (double **p, integer ni, integer i, integer k, double a)
  @          reads p[..][..];
  @ } */




/*@ 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);
  @ ensures \round_error(\result) <= 14*0x1.p-52 && \abs(\result) <= 2;
  @ */
double p_zero(double xs, double l, double x);


/*@ requires ni >= 2 && nk >= 2
  @          && l != 0
  @          && dx > 0. && dt > 0.  && v > 0.
  @          && \exact(dx) > 0. && \exact(dt) > 0.
  @          && \exact(v)==v
  @          && \abs(\exact(dx) - dx) / dx <= 0x1.p-53
  @          && \abs(\exact(dt) - dt) / dt <= 0x1.p-51
  @          && 3./5. <= \exact(dt)/\exact(dx) * \exact(v) <= 1-0x1.p-50
  @          && 0x1.p-1000 <= v <= 0x1.p1000
  @          && ni <= 0x1.p64 && nk <= 4194304
  @          && \exact(dx) <= 1;
  @
  @
  @ ensures \forall int i; \forall int k;
  @     0 <= i <= ni ==> 0 <= k <= nk ==> 
  @     \round_error(\result[i][k]) <= 78./2*0x1.p-52*(k+1)*(k+2); 
  @ */

double **forward_prop(int ni, int nk, double dx, double dt, double v,
                      double xs, double l) {

  /* output variable */
  double **p;
  
  /* local variables */
  int i, k;
  double a1, a, dp;

  
  /* assemble diagonal of stiffness matrix M[v]^{-1}A */
  a1 =  dt/dx*v;
  a  = a1*a1;
  
  /*@ assert 1./4 <= a <= 1 && 
    @         0 < \exact(a) <= 1 &&
    @         \round_error(a) <= 0x1.p-49;
    @  */

  p = array2d_alloc(ni+1, nk+1);
  /* initial conditions (includes boundary conditions) */

  p[0][0]=0.;
  /*@ loop invariant 1 <= i <= ni && analytic_error(p,ni,i-1,0,a);
    @ loop variant ni-i; */
  for (i=1; i<ni; i++) {
    p[i][0] = p_zero(xs, l, i*dx);
  }
  p[ni][0] =0.;

  /*@ assert analytic_error(p,ni,ni,0,a); */


  /* initial time derivative is supposed null + boundary conditions */

  p[0][1] = 0.;

  /*@ loop invariant 1 <= i <= ni && analytic_error(p,ni,i-1,1,a);
    @ loop variant ni-i; */
  for (i=1; i<ni; i++) {
    dp = p[i+1][0] - 2.*p[i][0] + p[i-1][0];
    p[i][1] = p[i][0] + 0.5*a*dp;
  }
  p[ni][1] = 0.;
  
  /*@ assert analytic_error(p,ni,ni,1,a); */


  /* propagation = time loop */
  /*@ loop invariant 1 <= k <= nk && analytic_error(p,ni,ni,k,a);
    @ loop variant nk-k; */
  for (k=1; k<nk; k++) {
    /* boundary condition 1 */
    p[0][k+1] = 0.;

    /* time iteration = space loop */
    /*@ loop invariant 1 <= i <= ni && analytic_error(p,ni,i-1,k+1,a);
      @ loop variant ni-i; */
 
    for (i=1; i<ni; i++) {
      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;
    }
    
    /* boundary condition 2 */
    p[ni][k+1] = 0.;
    /*@ assert analytic_error(p,ni,ni,k+1,a); */
    
  }

  return p;

}