#include <stdio.h>
#include <math.h>
#define EPSILON 1e-7
#define MAX_ITER 100
// 【確定版】プログラムが算出した正しい特性方程式 f(λ)
double f(double lambda) {
return pow(lambda
, 4) + pow(lambda
, 3) - 54.0 * pow(lambda
, 2) - 224.0 * lambda
- 345.0; }
// 【確定版】正しい導関数 f'(λ)
double df(double lambda) {
return 4.0 * pow(lambda
, 3) + 3.0 * pow(lambda
, 2) - 108.0 * lambda
- 224.0; }
void solve_newton(int attempt, double lambda0) {
double lambda = lambda0;
int iter = 0;
while (iter <= MAX_ITER) {
double f_val = f(lambda);
double df_val = df(lambda);
if (fabs(f_val
) < EPSILON
) { printf("試行 #%d (初期値 %5.1f) -> 収束成功 (%2d回) | 固有値解 λ = %.10f\n", attempt, lambda0, iter, lambda);
return;
}
if (fabs(df_val
) < 1e-12) { return; // 傾きが平らになった場合は終了
}
lambda -= f_val / df_val;
iter++;
}
}
int main() {
// 4つの解が重複せずにそれぞれ別の谷(根)に落ちるように調整した初期値
double initial_lambda[] = {10.0, -2.0, -3.2, -5.0};
printf("=== ニュートン・ラフソン法による実数解(固有値)の探索結果 ===\n"); for (int i = 0; i < 4; i++) {
solve_newton(i + 1, initial_lambda[i]);
}
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgojZGVmaW5lIEVQU0lMT04gMWUtNyAgIAojZGVmaW5lIE1BWF9JVEVSIDEwMAoKLy8g44CQ56K65a6a54mI44CR44OX44Ot44Kw44Op44Og44GM566X5Ye644GX44Gf5q2j44GX44GE54m55oCn5pa556iL5byPIGYozrspCmRvdWJsZSBmKGRvdWJsZSBsYW1iZGEpIHsKICAgIHJldHVybiBwb3cobGFtYmRhLCA0KSArIHBvdyhsYW1iZGEsIDMpIC0gNTQuMCAqIHBvdyhsYW1iZGEsIDIpIC0gMjI0LjAgKiBsYW1iZGEgLSAzNDUuMDsKfQoKLy8g44CQ56K65a6a54mI44CR5q2j44GX44GE5bCO6Zai5pWwIGYnKM67KQpkb3VibGUgZGYoZG91YmxlIGxhbWJkYSkgewogICAgcmV0dXJuIDQuMCAqIHBvdyhsYW1iZGEsIDMpICsgMy4wICogcG93KGxhbWJkYSwgMikgLSAxMDguMCAqIGxhbWJkYSAtIDIyNC4wOwp9Cgp2b2lkIHNvbHZlX25ld3RvbihpbnQgYXR0ZW1wdCwgZG91YmxlIGxhbWJkYTApIHsKICAgIGRvdWJsZSBsYW1iZGEgPSBsYW1iZGEwOwogICAgaW50IGl0ZXIgPSAwOwoKICAgIHdoaWxlIChpdGVyIDw9IE1BWF9JVEVSKSB7CiAgICAgICAgZG91YmxlIGZfdmFsID0gZihsYW1iZGEpOwogICAgICAgIGRvdWJsZSBkZl92YWwgPSBkZihsYW1iZGEpOwoKICAgICAgICBpZiAoZmFicyhmX3ZhbCkgPCBFUFNJTE9OKSB7CiAgICAgICAgICAgIHByaW50Zigi6Kmm6KGMICMlZCAo5Yid5pyf5YCkICU1LjFmKSAtPiDlj47mnZ/miJDlip8gKCUyZOWbnikgfCDlm7rmnInlgKTop6MgzrsgPSAlLjEwZlxuIiwgCiAgICAgICAgICAgICAgICAgICBhdHRlbXB0LCBsYW1iZGEwLCBpdGVyLCBsYW1iZGEpOwogICAgICAgICAgICByZXR1cm47CiAgICAgICAgfQoKICAgICAgICBpZiAoZmFicyhkZl92YWwpIDwgMWUtMTIpIHsKICAgICAgICAgICAgcmV0dXJuOyAvLyDlgr7jgY3jgYzlubPjgonjgavjgarjgaPjgZ/loLTlkIjjga/ntYLkuoYKICAgICAgICB9CgogICAgICAgIGxhbWJkYSAtPSBmX3ZhbCAvIGRmX3ZhbDsKICAgICAgICBpdGVyKys7CiAgICB9Cn0KCmludCBtYWluKCkgewogICAgLy8gNOOBpOOBruino+OBjOmHjeikh+OBm+OBmuOBq+OBneOCjOOBnuOCjOWIpeOBruiwt++8iOague+8ieOBq+iQveOBoeOCi+OCiOOBhuOBq+iqv+aVtOOBl+OBn+WIneacn+WApAogICAgZG91YmxlIGluaXRpYWxfbGFtYmRhW10gPSB7MTAuMCwgLTIuMCwgLTMuMiwgLTUuMH07CgogICAgcHJpbnRmKCI9PT0g44OL44Ol44O844OI44Oz44O744Op44OV44K944Oz5rOV44Gr44KI44KL5a6f5pWw6Kej77yI5Zu65pyJ5YCk77yJ44Gu5o6i57Si57WQ5p6cID09PVxuIik7CiAgICBmb3IgKGludCBpID0gMDsgaSA8IDQ7IGkrKykgewogICAgICAgIHNvbHZlX25ld3RvbihpICsgMSwgaW5pdGlhbF9sYW1iZGFbaV0pOwogICAgfQoKICAgIHJldHVybiAwOwp9Cg==