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