#include <stdio.h>
#include <math.h>
#define EPS 0.00000001
#define N 1000
typedef struct {
double ans;
int iter;
} Ans;
double f(double x)
{
return 0.5 - x
+ 0.2 * sin(x
); }
/* f'(x) */
double f2(double x)
{
return -1.0 + 0.2 * cos(x
); }
double g(double x)
{
return 0.5 + 0.2 * sin(x
); }
Ans linear_iteration_method(double x, double (*f)(double))
{
double x0;
int i = 0;
Ans ans;
do {
x0 = x;
x = f(x);
i++;
} while (fabs(x
- x0
) > EPS
);
ans.ans = x;
ans.iter = i;
return ans;
}
Ans aitken_deltasquare_process(double a[], double b[])
{
int i;
double dx1, dx2;
Ans ans;
for (i = 0; i < N - 2; i++) {
dx1 = a[i + 1] - a[i];
dx2 = a[i + 2] - a[i + 1];
b[i] = a[i + 2] - dx2 * dx2 / (dx2 - dx1);
if (fabs(b
[i
] - a
[i
+ 2]) <= EPS
) break;
}
ans.ans = b[i];
ans.iter = i + 1;
return ans;
}
Ans newton_method(double x, double (*f)(double), double (*f2)(double))
{
int i = 0;
double x0;
Ans ans;
do {
x0 = x;
x = x - f(x) / f2(x);
i++;
} while (fabs(x0
- x
) > EPS
);
ans.ans = x;
ans.iter = i;
return ans;
}
int main(void)
{
int i;
double x0 = 0.2;
Ans ans;
static double a[N], b[N];
ans = linear_iteration_method(x0, g);
printf("solution = %.8f, iteration = %d\n", ans.
ans, ans.
iter);
/* 元数列の生成 */
a[0] = x0;
for (i = 1; i < N; i++)
a[i] = g(a[i - 1]);
ans = aitken_deltasquare_process(a, b);
printf("solution = %.8f, iteration = %d\n", ans.
ans, ans.
iter);
ans = newton_method(x0, f, f2);
printf("solution = %.8f, iteration = %d\n", ans.
ans, ans.
iter);
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgojZGVmaW5lIEVQUyAwLjAwMDAwMDAxCiNkZWZpbmUgTiAgIDEwMDAKCnR5cGVkZWYgc3RydWN0IHsKICBkb3VibGUgYW5zOwogIGludCBpdGVyOwp9IEFuczsKCmRvdWJsZSBmKGRvdWJsZSB4KQp7CiAgcmV0dXJuIDAuNSAtIHggKyAwLjIgKiBzaW4oeCk7Cn0KCi8qIGYnKHgpICovCmRvdWJsZSBmMihkb3VibGUgeCkKewogIHJldHVybiAtMS4wICsgMC4yICogY29zKHgpOwp9Cgpkb3VibGUgZyhkb3VibGUgeCkKewogIHJldHVybiAwLjUgKyAwLjIgKiBzaW4oeCk7Cn0KCkFucyBsaW5lYXJfaXRlcmF0aW9uX21ldGhvZChkb3VibGUgeCwgZG91YmxlICgqZikoZG91YmxlKSkKewogIGRvdWJsZSB4MDsKICBpbnQgaSA9IDA7CiAgQW5zIGFuczsKICAKICBkbyB7CiAgICB4MCA9IHg7CiAgICB4ID0gZih4KTsKICAgIGkrKzsKICB9IHdoaWxlIChmYWJzKHggLSB4MCkgPiBFUFMpOwogIAogIGFucy5hbnMgPSB4OwogIGFucy5pdGVyID0gaTsKICAKICByZXR1cm4gYW5zOwp9CgpBbnMgYWl0a2VuX2RlbHRhc3F1YXJlX3Byb2Nlc3MoZG91YmxlIGFbXSwgZG91YmxlIGJbXSkKewogIGludCBpOwogIGRvdWJsZSBkeDEsIGR4MjsKICBBbnMgYW5zOwoKICBmb3IgKGkgPSAwOyBpIDwgTiAtIDI7IGkrKykgewogICAgZHgxID0gYVtpICsgMV0gLSBhW2ldOwogICAgZHgyID0gYVtpICsgMl0gLSBhW2kgKyAxXTsKICAgIGJbaV0gPSBhW2kgKyAyXSAtIGR4MiAqIGR4MiAvIChkeDIgLSBkeDEpOwogICAgaWYgKGZhYnMoYltpXSAtIGFbaSArIDJdKSA8PSBFUFMpCiAgICAgIGJyZWFrOwogIH0KCiAgYW5zLmFucyA9IGJbaV07CiAgYW5zLml0ZXIgPSBpICsgMTsKCiAgcmV0dXJuIGFuczsKfQoKQW5zIG5ld3Rvbl9tZXRob2QoZG91YmxlIHgsIGRvdWJsZSAoKmYpKGRvdWJsZSksIGRvdWJsZSAoKmYyKShkb3VibGUpKQp7CiAgaW50IGkgPSAwOwogIGRvdWJsZSB4MDsKICBBbnMgYW5zOwogIAogIGRvIHsKICAgIHgwID0geDsKICAgIHggPSB4IC0gZih4KSAvIGYyKHgpOwogICAgaSsrOwogIH0gd2hpbGUgKGZhYnMoeDAgLSB4KSA+IEVQUyk7CiAgCiAgYW5zLmFucyA9IHg7CiAgYW5zLml0ZXIgPSBpOwogIAogIHJldHVybiBhbnM7Cn0KCmludCBtYWluKHZvaWQpCnsKICBpbnQgaTsKICBkb3VibGUgeDAgPSAwLjI7CiAgQW5zIGFuczsKICBzdGF0aWMgZG91YmxlIGFbTl0sIGJbTl07CiAgCiAgcHV0cygi57ea5b2i5Y+N5b6p5rOVIik7CiAgYW5zID0gbGluZWFyX2l0ZXJhdGlvbl9tZXRob2QoeDAsIGcpOwogIHByaW50Zigic29sdXRpb24gPSAlLjhmLCBpdGVyYXRpb24gPSAlZFxuIiwgYW5zLmFucywgYW5zLml0ZXIpOwogIAogIHB1dHMoIuOCqOOCpOODiOOCseODs+OBruODh+ODq+OCvzLkuZfms5UiKTsKICAvKiDlhYPmlbDliJfjga7nlJ/miJAgKi8KICBhWzBdID0geDA7CiAgZm9yIChpID0gMTsgaSA8IE47IGkrKykKICAgIGFbaV0gPSBnKGFbaSAtIDFdKTsKICBhbnMgPSBhaXRrZW5fZGVsdGFzcXVhcmVfcHJvY2VzcyhhLCBiKTsKICBwcmludGYoInNvbHV0aW9uID0gJS44ZiwgaXRlcmF0aW9uID0gJWRcbiIsIGFucy5hbnMsIGFucy5pdGVyKTsKCiAgcHV0cygi44OL44Ol44O844OI44Oz5rOVIik7CiAgYW5zID0gbmV3dG9uX21ldGhvZCh4MCwgZiwgZjIpOwogIHByaW50Zigic29sdXRpb24gPSAlLjhmLCBpdGVyYXRpb24gPSAlZFxuIiwgYW5zLmFucywgYW5zLml0ZXIpOwoKICByZXR1cm4gMDsKfQo=