| 1 | /* Discrete derivative function using central differencing, which is
|
|---|
| 2 | * second-order accurate, except at the two end-points (where it is
|
|---|
| 3 | * first-order).
|
|---|
| 4 | * This version is incorrect, because it asserts third-order accuracy.
|
|---|
| 5 | * To attempt to verify with CIVL, type:
|
|---|
| 6 | * civl verify derivativeBad.cvl
|
|---|
| 7 | */
|
|---|
| 8 | $input double dx; // delta x
|
|---|
| 9 | $assume(0<dx && dx<1);
|
|---|
| 10 | $input int num_elements;
|
|---|
| 11 | $assume(num_elements >= 1);
|
|---|
| 12 | $input double in[num_elements];
|
|---|
| 13 | double out[num_elements];
|
|---|
| 14 | // the following says rho is a function from R to R which has 3 continuous
|
|---|
| 15 | // derivatives in the closed interval [-1,1]:
|
|---|
| 16 | $abstract $differentiable(3, [-1,1]) $real rho($real x);
|
|---|
| 17 |
|
|---|
| 18 | /* Computes discrete derivative by central differencing.
|
|---|
| 19 | * Right end-point is computed by backwards differencing.
|
|---|
| 20 | * Left end-point is computed by forward differencing.
|
|---|
| 21 | */
|
|---|
| 22 | void differentiate(int n, double y[], double h, double result[]) {
|
|---|
| 23 | $assume(n*h<=1);
|
|---|
| 24 | $assume($forall (int i : 0..n-1) y[i] == rho(i*h));
|
|---|
| 25 | /*@ loop invariant 1<=i && i<=n-1;
|
|---|
| 26 | @ loop invariant $forall (int j : 1..i-1) result[j] == (y[j+1]-y[j-1])/(2*h);
|
|---|
| 27 | @ loop assigns i, result[1..n-2];
|
|---|
| 28 | @*/
|
|---|
| 29 | for (int i=1; i<n-1; i++) {
|
|---|
| 30 | result[i] = (y[i+1]-y[i-1])/(2*h);
|
|---|
| 31 | }
|
|---|
| 32 | result[0] = (y[1]-y[0])/h;
|
|---|
| 33 | result[n-1] = (y[n-1] - y[n-2])/h;
|
|---|
| 34 | $assert($uniform (int i : 1..n-2) result[i]-$D[rho,{x,1}](i*h) == $O(h*h*h));
|
|---|
| 35 | }
|
|---|
| 36 |
|
|---|
| 37 | int main() {
|
|---|
| 38 | differentiate(num_elements, in, dx, out);
|
|---|
| 39 | }
|
|---|