Discretization of the 1D acoustic wave equation
This is a numerical analysis example that computes the discretization of the 1D acoustic wave equation.
Authors: Sylvie Boldo
Topics: Floating-Point Computations / Mathematics
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)
@ 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;
}