#include <stdio.h>
#include <math.h>
#define EPSILON 1e-7
#define MAX_ITER 100
// 1つ目の行列の特性方程式 f(λ)
double f(double lambda) {
return pow(lambda
, 4) + pow(lambda
, 3) - 88.0 * pow(lambda
, 2) - 132.0 * lambda
+ 603.0; }
// 導関数 f'(λ)
double df(double lambda) {
return 4.0 * pow(lambda
, 3) + 3.0 * pow(lambda
, 2) - 176.0 * lambda
- 132.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) { printf("試行 #%d (初期値 %5.1f) -> 導関数が0に近いため失敗\n", attempt
, lambda0
); return;
}
lambda -= f_val / df_val;
iter++;
}
printf("試行 #%d (初期値 %5.1f) -> %d回以内に収束せず失敗\n", attempt
, lambda0
, MAX_ITER
); }
int main() {
// 4つの異なる初期値からそれぞれの根(固有値)を直接探し出す
double initial_lambda[] = {8.0, 5.0, -1.0, -9.0};
printf("=== ニュートン・ラフソン法による実数解(固有値)の探索結果 ===\n"); for (int i = 0; i < 4; i++) {
solve_newton(i + 1, initial_lambda[i]);
}
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgojZGVmaW5lIEVQU0lMT04gMWUtNyAgIAojZGVmaW5lIE1BWF9JVEVSIDEwMAoKLy8gMeOBpOebruOBruihjOWIl+OBrueJueaAp+aWueeoi+W8jyBmKM67KQpkb3VibGUgZihkb3VibGUgbGFtYmRhKSB7CiAgICByZXR1cm4gcG93KGxhbWJkYSwgNCkgKyBwb3cobGFtYmRhLCAzKSAtIDg4LjAgKiBwb3cobGFtYmRhLCAyKSAtIDEzMi4wICogbGFtYmRhICsgNjAzLjA7Cn0KCi8vIOWwjumWouaVsCBmJyjOuykKZG91YmxlIGRmKGRvdWJsZSBsYW1iZGEpIHsKICAgIHJldHVybiA0LjAgKiBwb3cobGFtYmRhLCAzKSArIDMuMCAqIHBvdyhsYW1iZGEsIDIpIC0gMTc2LjAgKiBsYW1iZGEgLSAxMzIuMDsKfQoKdm9pZCBzb2x2ZV9uZXd0b24oaW50IGF0dGVtcHQsIGRvdWJsZSBsYW1iZGEwKSB7CiAgICBkb3VibGUgbGFtYmRhID0gbGFtYmRhMDsKICAgIGludCBpdGVyID0gMDsKCiAgICB3aGlsZSAoaXRlciA8PSBNQVhfSVRFUikgewogICAgICAgIGRvdWJsZSBmX3ZhbCA9IGYobGFtYmRhKTsKICAgICAgICBkb3VibGUgZGZfdmFsID0gZGYobGFtYmRhKTsKCiAgICAgICAgLy8g5Y+O5p2f5Yik5a6a44Gu44OB44Kn44OD44KvCiAgICAgICAgaWYgKGZhYnMoZl92YWwpIDwgRVBTSUxPTikgewogICAgICAgICAgICBwcmludGYoIuippuihjCAjJWQgKOWIneacn+WApCAlNS4xZikgLT4g5Y+O5p2f5oiQ5YqfICglMmTlm54pIHwg5Zu65pyJ5YCk6KejIM67ID0gJS4xMGZcbiIsIAogICAgICAgICAgICAgICAgICAgYXR0ZW1wdCwgbGFtYmRhMCwgaXRlciwgbGFtYmRhKTsKICAgICAgICAgICAgcmV0dXJuOwogICAgICAgIH0KCiAgICAgICAgaWYgKGZhYnMoZGZfdmFsKSA8IDFlLTEyKSB7CiAgICAgICAgICAgIHByaW50Zigi6Kmm6KGMICMlZCAo5Yid5pyf5YCkICU1LjFmKSAtPiDlsI7plqLmlbDjgYww44Gr6L+R44GE44Gf44KB5aSx5pWXXG4iLCBhdHRlbXB0LCBsYW1iZGEwKTsKICAgICAgICAgICAgcmV0dXJuOwogICAgICAgIH0KCiAgICAgICAgbGFtYmRhIC09IGZfdmFsIC8gZGZfdmFsOwogICAgICAgIGl0ZXIrKzsKICAgIH0KICAgIHByaW50Zigi6Kmm6KGMICMlZCAo5Yid5pyf5YCkICU1LjFmKSAtPiAlZOWbnuS7peWGheOBq+WPjuadn+OBm+OBmuWkseaVl1xuIiwgYXR0ZW1wdCwgbGFtYmRhMCwgTUFYX0lURVIpOwp9CgppbnQgbWFpbigpIHsKICAgIC8vIDTjgaTjga7nlbDjgarjgovliJ3mnJ/lgKTjgYvjgonjgZ3jgozjgZ7jgozjga7moLnvvIjlm7rmnInlgKTvvInjgpLnm7TmjqXmjqLjgZflh7rjgZkKICAgIGRvdWJsZSBpbml0aWFsX2xhbWJkYVtdID0gezguMCwgNS4wLCAtMS4wLCAtOS4wfTsKCiAgICBwcmludGYoIj09PSDjg4vjg6Xjg7zjg4jjg7Pjg7vjg6njg5Xjgr3jg7Pms5Xjgavjgojjgovlrp/mlbDop6PvvIjlm7rmnInlgKTvvInjga7mjqLntKLntZDmnpwgPT09XG4iKTsKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgNDsgaSsrKSB7CiAgICAgICAgc29sdmVfbmV3dG9uKGkgKyAxLCBpbml0aWFsX2xhbWJkYVtpXSk7CiAgICB9CgogICAgcmV0dXJuIDA7Cn0=