1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
| #include <math.h>
#include <stdio.h>
#define N 10
void richardson(double* h_k, double* a_k, double a, double* r_h_k, double* r_a_k) {
int i;
int N_a_k = 0;
double ratio, demo_D, a_k_richardson, p, c, distances[N];
double r[N*2];
/* avoid richarson extrapolation breaks down, truncate
for the first k where a_{k} = a_{k + 1}. We compute
the usable length of a_k */
for (i = 1; i < N; i++) {
if (fabs(a_k[i] - a_k[i-1]) < 1*pow(10,-14)) {
break;
}
N_a_k++;
}
ratio = h_k[1] / h_k[0];
demo_D = h_k[0] / ratio;
printf("Estimate order p and c_p such that a = a_k + c_p pow(h_k,p) + O(pow(h_k,[p+1]))\n--with the use of 3 subsequent a_k s:\n");
printf("Input is sequence a_k and h_k = %f * %f^k.\n", demo_D, ratio);
printf("+----------------+--------+--------+-----------------++-----------------+\n");
printf("| a_k | p | c_p | a_k_richardson || a |\n");
printf("+----------------+--------+--------+-----------------++-----------------+\n");
for (i = 0; i < N_a_k; i++){
if (i < 2) {
printf("|%15.11f | | | || |\n",a_k[i]);}
else
{ p = log((a_k[i] - a_k[i - 1]) / (a_k[i - 1] - a_k[i - 2])) / log(ratio);
c = (a_k[i] - a_k[i - 1]) / (pow(h_k[i - 1],p) - pow(h_k[i],p));
a_k_richardson = a_k[i] + c * pow(h_k[i],p);
printf("|%15.11f | %6.3f | %6.3f | %15.11f || %15.11f |\n", a_k[i], p, c, a_k_richardson, a - a_k[i]);
}
}
for (i = 0; i < N*2; i++)
r[i] = 0;
for (i = 0; i < N_a_k; i++) {
if (i < 1) {
printf("|%15.11f | | || |\n",a_k[i]); }
else {
c = (a_k[i] - a_k[i - 1])/(pow(h_k[i - 1],p) - pow(h_k[i],p));
r[i*2] = r[(i-1)*2 + 1];
r[i*2 + 1] = (a_k[i - 1] + c*pow(h_k[i - 1],p))/2 + (a_k[i] + c*pow(h_k[i],p))/2;
printf("|%15.11f | %15.11f | %6.3f || %15.11f |\n", a_k[i], r[i - 1], c, a - r[i - 1]);
}
}
printf("+----------------+-----------------+--------++-----------------+\n");
for (i = 0; i < N*2; i++)
{
r_a_k[i] = r[i];
r_h_k[i] = h_k[i];
}
printf("+----------------+--------+--------+-----------------++-----------------+\n");
for (i = 0; i < N_a_k; i++)
distances[i] = fabs(a_k[i] - a_k[N_a_k-1]);
p = round(p);
printf("The convergence order is p = %f, i.e.: a = a_k + a_%f h_k^%f + O(h_k^%f).\n", p, p, p, p + 1);
printf("Richardson extrapolants pow(a_k,1) with order %d:\n", p);
printf("+----------------+-----------------+--------++-----------------+\n");
printf("| a_k | a_k^(1) | c_p || a |\n");
printf("+----------------+-----------------+--------++-----------------+\n");
}
int main(void) {
double h_k[N], a_k[N], r_h_k[N*2], r_a_k[N*2];
h_k[0] = 1.0;
h_k[1] = 0.5;
/* ... fill input arrays */
richardson(h_k, a_k, 2.0, r_h_k, r_a_k);
system("pause");
return 0;
} |
