| [15d19f3] | 1 | /* A discrete approximation to the second derivative of a function from R to R.
|
|---|
| 2 | * It is second-order accurate, except at the two end-points.
|
|---|
| 3 | * This version is incorrect because it claims third-order accuracy.
|
|---|
| 4 | * To attempt to verify with CIVL, type:
|
|---|
| 5 | * civl verify secondDerivativeBad.cvl
|
|---|
| [965b371] | 6 | *
|
|---|
| 7 | * Note: based on Quarteroni, Sacco, Saleri. "Numerical Mathematics" 2nd ed. sec 10.10.1
|
|---|
| [15d19f3] | 8 | */
|
|---|
| [965b371] | 9 |
|
|---|
| [15d19f3] | 10 | $input double dx;
|
|---|
| 11 | $assume(0<dx && dx<1);
|
|---|
| 12 | $input int num_elements;
|
|---|
| 13 | $assume(num_elements > 0);
|
|---|
| 14 | $input double in[num_elements];
|
|---|
| 15 | double out[num_elements];
|
|---|
| 16 | // assume rho:R->R has 4 continuous derivatives on [-1,1]:
|
|---|
| 17 | $abstract $differentiable(4, [-1,1]) $real rho($real x);
|
|---|
| [965b371] | 18 |
|
|---|
| [15d19f3] | 19 | void secondDerivative(int n, double y[], double h, double result[]) {
|
|---|
| 20 | $assume($forall (int i:0..n-1) y[i] == rho(i*h));
|
|---|
| 21 | /*@ loop invariant 1<=i && i<=n-1;
|
|---|
| 22 | @ loop invariant $forall (int j:1..i-1)
|
|---|
| 23 | @ result[j] == (y[j+1] - 2*y[j] + y[j-1])/(h*h);
|
|---|
| 24 | @ loop assigns i, result[1..n-2];
|
|---|
| 25 | @*/
|
|---|
| 26 | for (int i = 1; i < n-1; i++) {
|
|---|
| 27 | result[i] = (y[i+1] - 2*y[i] + y[i-1])/(h*h);
|
|---|
| 28 | }
|
|---|
| 29 | result[0] = (y[2] - 2*y[1] + y[0])/h;
|
|---|
| 30 | result[n-1] = (y[n-3] - 2*y[n-2] - y[n-1])/h;
|
|---|
| 31 | $assert($uniform (int i:1..n-2) result[i] - $D[rho,{x,2}](i*h) == $O(h*h*h));
|
|---|
| [965b371] | 32 | }
|
|---|
| 33 |
|
|---|
| 34 | void main() {
|
|---|
| [15d19f3] | 35 | secondDerivative(num_elements, in, dx, out);
|
|---|
| [965b371] | 36 | }
|
|---|