#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;
printf("\n--- 試行 #%d: (初期値 lambda0 = %.2f) ---\n", attempt
, lambda0
); printf("%-5s %-15s\n", "反復", "lambda"); printf("-------------------------\n");
while (iter <= MAX_ITER) {
double f_val = f(lambda);
double df_val = df(lambda);
// lambda の値のみを綺麗に出力
printf("%-5d %-15.10f\n", iter
, lambda
);
// 収束判定
if (fabs(f_val
) < EPSILON
) { printf(">> 収束しました(反復回数: %d回)\n", iter
); printf(">> 得られた固有値 (lambda): %.10f\n", lambda
); return;
}
// 導関数がゼロに近すぎる場合の安全対策
if (fabs(df_val
) < 1e-12) { printf(">> 失敗: 導関数が0に近づきすぎました。\n"); return;
}
// ニュートン・ラフソン法の更新式: λ = λ - f(λ)/f'(λ)
lambda -= f_val / df_val;
iter++;
}
printf(">> 失敗: %d回以内に収束しませんでした。\n", MAX_ITER
); }
int main() {
// すべての実数解(固有値)を見つけられるように複数の初期値を設定
// 8.0 からスタートすると、べき乗法で求めた最大固有値に収束します
double initial_lambda[] = {8.0, 5.0, -1.0, -9.0};
for (int i = 0; i < 4; i++) {
solve_newton(i + 1, initial_lambda[i]);
}
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgojZGVmaW5lIEVQU0lMT04gMWUtNyAgIC8vIOOBueOBjeS5l+azleOBqOWQiOOCj+OBm+OBn+WPjuadn+WIpOWumuadoeS7tgojZGVmaW5lIE1BWF9JVEVSIDEwMAoKLy8gMeOBpOebruOBruihjOWIl+OBrueJueaAp+aWueeoi+W8jyBmKM67KQpkb3VibGUgZihkb3VibGUgbGFtYmRhKSB7CiAgICByZXR1cm4gcG93KGxhbWJkYSwgNCkgKyBwb3cobGFtYmRhLCAzKSAtIDg4LjAgKiBwb3cobGFtYmRhLCAyKSAtIDEzMi4wICogbGFtYmRhICsgNjAzLjA7Cn0KCi8vIOWwjumWouaVsCBmJyjOuykKZG91YmxlIGRmKGRvdWJsZSBsYW1iZGEpIHsKICAgIHJldHVybiA0LjAgKiBwb3cobGFtYmRhLCAzKSArIDMuMCAqIHBvdyhsYW1iZGEsIDIpIC0gMTc2LjAgKiBsYW1iZGEgLSAxMzIuMDsKfQoKdm9pZCBzb2x2ZV9uZXd0b24oaW50IGF0dGVtcHQsIGRvdWJsZSBsYW1iZGEwKSB7CiAgICBkb3VibGUgbGFtYmRhID0gbGFtYmRhMDsKICAgIGludCBpdGVyID0gMDsKCiAgICBwcmludGYoIlxuLS0tIOippuihjCAjJWQ6ICjliJ3mnJ/lgKQgbGFtYmRhMCA9ICUuMmYpIC0tLVxuIiwgYXR0ZW1wdCwgbGFtYmRhMCk7CiAgICBwcmludGYoIiUtNXMgJS0xNXNcbiIsICLlj43lvqkiLCAibGFtYmRhIik7CiAgICBwcmludGYoIi0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS1cbiIpOwoKICAgIHdoaWxlIChpdGVyIDw9IE1BWF9JVEVSKSB7CiAgICAgICAgZG91YmxlIGZfdmFsID0gZihsYW1iZGEpOwogICAgICAgIGRvdWJsZSBkZl92YWwgPSBkZihsYW1iZGEpOwoKICAgICAgICAvLyBsYW1iZGEg44Gu5YCk44Gu44G/44KS57a66bqX44Gr5Ye65YqbCiAgICAgICAgcHJpbnRmKCIlLTVkICUtMTUuMTBmXG4iLCBpdGVyLCBsYW1iZGEpOwoKICAgICAgICAvLyDlj47mnZ/liKTlrpoKICAgICAgICBpZiAoZmFicyhmX3ZhbCkgPCBFUFNJTE9OKSB7CiAgICAgICAgICAgIHByaW50ZigiPj4g5Y+O5p2f44GX44G+44GX44Gf77yI5Y+N5b6p5Zue5pWwOiAlZOWbnu+8iVxuIiwgaXRlcik7CiAgICAgICAgICAgIHByaW50ZigiPj4g5b6X44KJ44KM44Gf5Zu65pyJ5YCkIChsYW1iZGEpOiAlLjEwZlxuIiwgbGFtYmRhKTsKICAgICAgICAgICAgcmV0dXJuOwogICAgICAgIH0KCiAgICAgICAgLy8g5bCO6Zai5pWw44GM44K844Ot44Gr6L+R44GZ44GO44KL5aC05ZCI44Gu5a6J5YWo5a++562WCiAgICAgICAgaWYgKGZhYnMoZGZfdmFsKSA8IDFlLTEyKSB7CiAgICAgICAgICAgIHByaW50ZigiPj4g5aSx5pWXOiDlsI7plqLmlbDjgYww44Gr6L+R44Gl44GN44GZ44GO44G+44GX44Gf44CCXG4iKTsKICAgICAgICAgICAgcmV0dXJuOwogICAgICAgIH0KCiAgICAgICAgLy8g44OL44Ol44O844OI44Oz44O744Op44OV44K944Oz5rOV44Gu5pu05paw5byPOiDOuyA9IM67IC0gZijOuykvZicozrspCiAgICAgICAgbGFtYmRhIC09IGZfdmFsIC8gZGZfdmFsOwogICAgICAgIGl0ZXIrKzsKICAgIH0KICAgIHByaW50ZigiPj4g5aSx5pWXOiAlZOWbnuS7peWGheOBq+WPjuadn+OBl+OBvuOBm+OCk+OBp+OBl+OBn+OAglxuIiwgTUFYX0lURVIpOwp9CgppbnQgbWFpbigpIHsKICAgIC8vIOOBmeOBueOBpuOBruWun+aVsOino++8iOWbuuacieWApO+8ieOCkuimi+OBpOOBkeOCieOCjOOCi+OCiOOBhuOBq+ikh+aVsOOBruWIneacn+WApOOCkuioreWumgogICAgLy8gOC4wIOOBi+OCieOCueOCv+ODvOODiOOBmeOCi+OBqOOAgeOBueOBjeS5l+azleOBp+axguOCgeOBn+acgOWkp+WbuuacieWApOOBq+WPjuadn+OBl+OBvuOBmQogICAgZG91YmxlIGluaXRpYWxfbGFtYmRhW10gPSB7OC4wLCA1LjAsIC0xLjAsIC05LjB9OwoKICAgIGZvciAoaW50IGkgPSAwOyBpIDwgNDsgaSsrKSB7CiAgICAgICAgc29sdmVfbmV3dG9uKGkgKyAxLCBpbml0aWFsX2xhbWJkYVtpXSk7CiAgICB9CgogICAgcmV0dXJuIDA7Cn0=