C言語のプログラムで質問です。
ヤコビ法で固有値・固有ベクトル絵を求めたいのですが、固有値は上手く出るのですが、
固有ベクトルがめちゃくちゃな値が出ます。何が間違っているのでしょうか教えてください。下のは一応自分で作ってみたプログラムです。
#include <stdio.h>
#include <math.h>
#define N 4
#define PI 3.14159
void jacobi(double a[N][N], double λ[N], double v[N][N]);
int main(void)
{
int i, j;
double a[N][N],
λ[N], v[N][N], x;
for(i=0;i<N;i++){
printf("\n a[%d][0] a[%d][1] a[%d][2] a[%d][3]=",i,i,i,i);
scanf("%lf %lf %lf %lf",&a[i][0],&a[i][1],&a[i][2],&a[i][3]);
}
jacobi(a, λ, v);
for(j=0; j<N; j++){
printf("\n固有値λ[%d]=%lf", j, λ[j]);
printf("\n固有ベクトル\n");
for(i=0; i<N; i++){
printf("%lf\n", v[i][j]);
}
}
return 0;
}
/* ヤコビ法による固有値計算 */
void jacobi(double a[N][N], double λ[N], double v[N][N])
{
int i, j, kmax=100, repeat, p, q;
double eps, c, s,theta, gmax;
double apq, app, aqq, apqmax, apj, aqj, vip, viq;
gmax=0.0;
for(i=0; i<N; i++){
s=0.0;
for(j=i+1; j<N; j++){
s += fabs(a[i][j]);
}
if(s>gmax) gmax=s;
}
eps=0.000001*gmax;
for(i=0; i<N; i++){
for(j=0; j<N; j++){
v[i][j]=0.0;
}
v[i][i]=1.0;
}
for(repeat=1; repeat<kmax; repeat++){
/* 収束判定 */
apqmax = 0.0;
for(p=0; p<N; p++){
for(q=0; q<N; q++){
if(p!=q){
apq=fabs(a[p][q]);
if(apq>apqmax) apqmax=apq;
}
}
}
if(apqmax<eps) break;
for(p=0; p<N-1; p++){
for(q=p+1; q<N; q++){
apq=a[p][q];
app=a[p][p];
aqq=a[q][q];
if(fabs(apqmax)<eps) break;
/*回転角計算*/
if(fabs(app-aqq)>=1.0e-15){
theta= 0.5*atan(2.0*apq/(app-aqq));
}else{
theta= PI/4.0;
}
c = cos(theta);
s = sin(theta);
a[p][p] = app*c*c + 2.0*apq*c*s + aqq*s*s;
a[q][q] = app*s*s - 2.0*apq*c*s + aqq*c*c;
a[p][q] = 0.0;
a[q][p] = 0.0;
for(j=0; j<N; j++){
if(j!=p && j!=q){
apj = a[p][j];
aqj = a[q][j];
a[p][j] = apj*c + aqj*s;
a[q][j] = -apj*s + aqj*c;
a[j][p] = a[p][j];
a[j][q] = a[q][j];
}
}
/*固有ベクトル*/
for(i=0; i<N; i++){
v[i][p] = v[i][p]*c+v[i][q]*s;
v[i][q] = -v[i][p]*s+v[i][q]*c;
}
}
}
eps=eps*1.05;
}
for(i=0; i<N; i++) λ[i] = a[i][i];
}
No.1
- 回答日時:
なんか何をやっているのかよくわからんプログラムになってるなぁ.... apqmax と eps の比較はなんでこんなに必要なんだろう.
参考URL:http://hooktail.org/computer/index.php?Jacobi%CB …
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- C言語・C++・C# 10個の実数に対する降順ソート結果を出力するプログラムを作りたいのですが、以下のプログラムをどう直せ 1 2022/07/09 22:16
- C言語・C++・C# C 言語の Gauss Jordan 法について 2 2022/12/28 11:16
- C言語・C++・C# c言語でユーザ関数を利用して複素数のべき乗と絶対値の数列を計算するプログラムが作りたいです。 3 2023/01/29 22:13
- C言語・C++・C# LU分解法のピボッティングについて(C言語/gcc-9) 3 2022/07/11 23:10
- C言語・C++・C# LU分解法のピボット選択機能実装について(C言語・gcc-9) 1 2022/07/22 15:20
- C言語・C++・C# Cのdoubleの浮動小数点表示について 3 2023/04/17 13:14
- C言語・C++・C# 並列プログラミングのπ計算について 1 2022/07/16 22:30
- C言語・C++・C# C言語のエラーについて 2 2022/07/11 13:56
- C言語・C++・C# プログラミングの授業の課題です 1 2023/01/17 22:15
- C言語・C++・C# C言語 プログラミング 4 2022/05/22 11:53
関連するカテゴリからQ&Aを探す
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
プログラムでの数字につく”f”の...
-
float型とdouble型の変数の違い...
-
doubleの変数にintとintの割り...
-
C言語 関数プロトタイプ宣言の...
-
C言語で台形公式を使った二重積...
-
関数におけるif文とreturn文に...
-
C言語(プログラミング)関連の質...
-
C言語を実行すると-infが出てき...
-
C 開放してるのにエラー(doubl...
-
数値を指数部と仮数部に分離したい
-
マチンの公式による円周率のプ...
-
至急です! マクロ定義で #defi...
-
C# 分秒表示ついて
-
c言語で、繰り返し文の中で、0....
-
C言語でdouble型の小数点の引き...
-
学校の課題で2次方程式のプログ...
-
ニュートン法
-
C言語で表記についの質問です
-
c言語の問題
-
2分法で方程式の複数の解を自...
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
プログラムでの数字につく”f”の...
-
float型とdouble型の変数の違い...
-
doubleの変数にintとintの割り...
-
C言語を実行すると-infが出てき...
-
C 開放してるのにエラー(doubl...
-
至急です! マクロ定義で #defi...
-
c言語で、繰り返し文の中で、0....
-
関数におけるif文とreturn文に...
-
C言語 関数プロトタイプ宣言の...
-
C言語初心者 構造体 課題について
-
C言語の型による処理速度の違い
-
Cで3乗根を求める方法
-
C言語で-23乗を取り扱うには
-
2分法で方程式の複数の解を自...
-
doubleは常に%lfとするべきなのか
-
c言語のコンパイルエラー canno...
-
C言語で直角三角形の斜辺を求め...
-
C言語のプログラムで#include<m...
-
int とdoubleの比較
-
C++で外積
おすすめ情報