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; }