## Discretization of the 1D acoustic wave equation

This is a numerical analysis example that computes the discretization of the 1D acoustic wave equation.

Auteurs: Sylvie Boldo

Catégories: Floating-Point Computations / Mathematics

Outils: Frama-C / Jessie / Coq

see also the index (by topic, by tool, by reference, by year)

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.
/*@ 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)
@
@  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 ==>
@
@ logic real sqr_norm_dx_conv_err{L}
@   (double **p, real dx, real dt, integer ni, integer i, integer k)
@ 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;

}