| 情報学部 | 菅沼ホーム | 目次 | 索引 |
/********************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* coded by Y.Suganuma */
/********************************************/
#include <stdio.h>
int Bairstow(int, int, double, double, double,
double *, double *, double *, double *, double *);
int Frame(int, int, double, double, double, double *, double *, double *,
double *, double *, double **, double **, double **);
int main()
{
double **A, **H1, **H2, *a, *b, *c, *rl, *im, p0, q0, eps;
int i1, ind, ct, n;
// データの設定
ct = 1000;
eps = 1.0e-10;
p0 = 0.0;
q0 = 0.0;
n = 3;
a = new double [n+1];
b = new double [n+1];
c = new double [n+1];
rl = new double [n];
im = new double [n];
A = new double * [n];
H1 = new double * [n];
H2 = new double * [n];
for (i1 = 0; i1 < n; i1++) {
A[i1] = new double [n];
H1[i1] = new double [n];
H2[i1] = new double [n];
}
A[0][0] = 7.0;
A[0][1] = 2.0;
A[0][2] = 1.0;
A[1][0] = 2.0;
A[1][1] = 1.0;
A[1][2] = -4.0;
A[2][0] = 1.0;
A[2][1] = -4.0;
A[2][2] = 2.0;
// 計算
ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2);
// 出力
if (ind > 0)
printf("収束しませんでした！\n");
else {
for (i1 = 0; i1 < n; i1++)
printf(" %f i %f\n", rl[i1], im[i1]);
}
delete [] a;
delete [] b;
delete [] c;
delete [] rl;
delete [] im;
for (i1 = 0; i1 < n; i1++) {
delete [] A[i1];
delete [] H1[i1];
delete [] H2[i1];
}
delete [] A;
delete [] H1;
delete [] H2;
return 0;
}
/*************************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b, c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* A : 行列 */
/* H1, H2 : 作業域（nxnの行列） */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
int Frame(int n, int ct, double eps, double p0, double q0, double *a, double *b,
double *c, double *rl, double *im, double **A, double **H1, double **H2)
{
int i1, i2, i3, i4, ind;
a[0] = 1.0;
// a1の計算
a[1] = 0.0;
for (i1 = 0; i1 < n; i1++)
a[1] -= A[i1][i1];
// a2の計算
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
H1[i1][i2] = A[i1][i2];
H1[i1][i1] += a[1];
}
a[2] = 0.0;
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
a[2] -= A[i1][i2] * H1[i2][i1];
}
a[2] *= 0.5;
// a3からanの計算
for (i1 = 3; i1 <= n; i1++) {
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++) {
H2[i2][i3] = 0.0;
for (i4 = 0; i4 < n; i4++)
H2[i2][i3] += A[i2][i4] * H1[i4][i3];
}
H2[i2][i2] += a[i1-1];
}
a[i1] = 0.0;
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++)
a[i1] -= A[i2][i3] * H2[i3][i2];
}
a[i1] /= i1;
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++)
H1[i2][i3] = H2[i2][i3];
}
}
// ベアストウ法
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im);
return ind;
}
/*************************************************/
/* 実係数代数方程式の解（ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b, c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
#include <math.h>
int Bairstow(int n, int ct, double eps, double p0, double q0,
double *a, double *b, double *c, double *rl, double *im)
{
double D, dp, dq, p1 = p0, p2 = 0.0, q1 = q0, q2 = 0.0;
int i1, ind = 0, count = 0;
/*
１次の場合
*/
if (n == 1) {
if (fabs(a[0]) < eps)
ind = 1;
else {
rl[0] = -a[1] / a[0];
im[0] = 0.0;
}
}
/*
２次の場合
*/
else if (n == 2) {
// １次式
if (fabs(a[0]) < eps) {
if (fabs(a[1]) < eps)
ind = 1;
else {
rl[0] = -a[2] / a[1];
im[0] = 0.0;
}
}
// ２次式
else {
D = a[1] * a[1] - 4.0 * a[0] * a[2];
if (D < 0.0) { // 虚数
D = sqrt(-D);
a[0] *= 2.0;
rl[0] = -a[1] / a[0];
rl[1] = -a[1] / a[0];
im[0] = D / a[0];
im[1] = -im[0];
}
else { // 実数
D = sqrt(D);
a[0] = 1.0 / (2.0 * a[0]);
rl[0] = a[0] * (-a[1] + D);
rl[1] = a[0] * (-a[1] - D);
im[0] = 0.0;
im[1] = 0.0;
}
}
}
// ３次以上の場合
else {
// 因数分解
ind = 1;
while (ind > 0 && count <= ct) {
for (i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
b[i1] = a[i1];
else if (i1 == 1)
b[i1] = a[i1] - p1 * b[i1-1];
else
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2];
}
for (i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
c[i1] = b[i1];
else if (i1 == 1)
c[i1] = b[i1] - p1 * c[i1-1];
else
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2];
}
D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1]);
if (fabs(D) < eps)
return ind;
else {
dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / D;
dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / D;
p2 = p1 + dp;
q2 = q1 + dq;
if (fabs(dp) < eps && fabs(dq) < eps)
ind = 0;
else {
count++;
p1 = p2;
q1 = q2;
}
}
}
if (ind == 0) {
// ２次方程式を解く
D = p2 * p2 - 4.0 * q2;
if (D < 0.0) { // 虚数
D = sqrt(-D);
rl[0] = -0.5 * p2;
rl[1] = -0.5 * p2;
im[0] = 0.5 * D;
im[1] = -im[0];
}
else { // 実数
D = sqrt(D);
rl[0] = 0.5 * (-p2 + D);
rl[1] = 0.5 * (-p2 - D);
im[0] = 0.0;
im[1] = 0.0;
}
// 残りの方程式を解く
n -= 2;
for (i1 = 0; i1 <= n; i1++)
a[i1] = b[i1];
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, &rl[2], &im[2]);
}
}
return ind;
}
/********************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* coded by Y.Suganuma */
/********************************************/
import java.io.*;
public class Test {
public static void main(String args[]) throws IOException
{
double A[][], H1[][], H2[][], a[], b[], c[], rl[], im[], p0, q0, eps;
int i1, ind, ct, n;
// データの設定
ct = 1000;
eps = 1.0e-10;
p0 = 0.0;
q0 = 0.0;
n = 3;
a = new double [n+1];
b = new double [n+1];
c = new double [n+1];
rl = new double [n];
im = new double [n];
A = new double [n][n];
H1 = new double [n][n];
H2 = new double [n][n];
A[0][0] = 7.0;
A[0][1] = 2.0;
A[0][2] = 1.0;
A[1][0] = 2.0;
A[1][1] = 1.0;
A[1][2] = -4.0;
A[2][0] = 1.0;
A[2][1] = -4.0;
A[2][2] = 2.0;
// 計算
ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2);
// 出力
if (ind > 0)
System.out.println("収束しませんでした！");
else {
for (i1 = 0; i1 < n; i1++)
System.out.println(" " + rl[i1] + " i " + im[i1]);
}
}
/*************************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b, c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* A : 行列 */
/* H1, H2 : 作業域（nxnの行列） */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
static int Frame(int n, int ct, double eps, double p0, double q0, double a[], double b[],
double c[], double rl[], double im[], double A[][], double H1[][], double H2[][])
{
int i1, i2, i3, i4, ind;
a[0] = 1.0;
// a1の計算
a[1] = 0.0;
for (i1 = 0; i1 < n; i1++)
a[1] -= A[i1][i1];
// a2の計算
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
H1[i1][i2] = A[i1][i2];
H1[i1][i1] += a[1];
}
a[2] = 0.0;
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
a[2] -= A[i1][i2] * H1[i2][i1];
}
a[2] *= 0.5;
// a3からanの計算
for (i1 = 3; i1 <= n; i1++) {
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++) {
H2[i2][i3] = 0.0;
for (i4 = 0; i4 < n; i4++)
H2[i2][i3] += A[i2][i4] * H1[i4][i3];
}
H2[i2][i2] += a[i1-1];
}
a[i1] = 0.0;
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++)
a[i1] -= A[i2][i3] * H2[i3][i2];
}
a[i1] /= i1;
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++)
H1[i2][i3] = H2[i2][i3];
}
}
// ベアストウ法
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);
return ind;
}
/*************************************************/
/* 実係数代数方程式の解（ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b,c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* p : 答えの位置 */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
static int Bairstow(int n, int ct, double eps, double p0, double q0,
double a[], double b[], double c[], double rl[], double im[], int p)
{
double D, dp, dq, p1 = p0, p2 = 0.0, q1 = q0, q2 = 0.0;
int i1, ind = 0, count = 0;
/*
１次の場合
*/
if (n == 1) {
if (Math.abs(a[0]) < eps)
ind = 1;
else {
rl[p] = -a[1] / a[0];
im[p] = 0.0;
}
}
/*
２次の場合
*/
else if (n == 2) {
// １次式
if (Math.abs(a[0]) < eps) {
if (Math.abs(a[1]) < eps)
ind = 1;
else {
rl[p] = -a[2] / a[1];
im[p] = 0.0;
}
}
// ２次式
else {
D = a[1] * a[1] - 4.0 * a[0] * a[2];
if (D < 0.0) { // 虚数
D = Math.sqrt(-D);
a[0] *= 2.0;
rl[p] = -a[1] / a[0];
rl[p+1] = -a[1] / a[0];
im[p] = D / a[0];
im[p+1] = -im[p];
}
else { // 実数
D = Math.sqrt(D);
a[0] = 1.0 / (2.0 * a[0]);
rl[p] = a[0] * (-a[1] + D);
rl[p+1] = a[0] * (-a[1] - D);
im[p] = 0.0;
im[p+1] = 0.0;
}
}
}
// ３次以上の場合
else {
// 因数分解
ind = 1;
while (ind > 0 && count <= ct) {
for (i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
b[i1] = a[i1];
else if (i1 == 1)
b[i1] = a[i1] - p1 * b[i1-1];
else
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2];
}
for (i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
c[i1] = b[i1];
else if (i1 == 1)
c[i1] = b[i1] - p1 * c[i1-1];
else
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2];
}
D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1]);
if (Math.abs(D) < eps)
return ind;
else {
dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / D;
dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / D;
p2 = p1 + dp;
q2 = q1 + dq;
if (Math.abs(dp) < eps && Math.abs(dq) < eps)
ind = 0;
else {
count++;
p1 = p2;
q1 = q2;
}
}
}
if (ind == 0) {
// ２次方程式を解く
D = p2 * p2 - 4.0 * q2;
if (D < 0.0) { // 虚数
D = Math.sqrt(-D);
rl[p] = -0.5 * p2;
rl[p+1] = -0.5 * p2;
im[p] = 0.5 * D;
im[p+1] = -im[p];
}
else { // 実数
D = Math.sqrt(D);
rl[p] = 0.5 * (-p2 + D);
rl[p+1] = 0.5 * (-p2 - D);
im[p] = 0.0;
im[p+1] = 0.0;
}
// 残りの方程式を解く
n -= 2;
for (i1 = 0; i1 <= n; i1++)
a[i1] = b[i1];
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p+2);
}
}
return ind;
}
}
<!DOCTYPE HTML>
<HTML>
<HEAD>
<TITLE>行列の固有値（フレーム法＋ベアストウ法）</TITLE>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
<SCRIPT TYPE="text/javascript">
function main()
{
// データの設定
let ct = parseInt(document.getElementById("trial").value);
let eps = 1.0e-10;
let p0 = parseFloat(document.getElementById("p0").value);
let q0 = parseFloat(document.getElementById("q0").value);
let n = parseInt(document.getElementById("order").value);
let a = new Array();
let b = new Array();
let c = new Array();
let rl = new Array();
let im = Array();
let A = new Array();
for (let i1 = 0; i1 < n; i1++)
A[i1] = new Array();
let s = (document.getElementById("ar").value).split(/ {1,}|\n{1,}/);
let k = 0;
for (let i1 = 0; i1 < n; i1++) {
for (let i2 = 0; i2 < n; i2++) {
A[i1][i2] = parseFloat(s[k]);
k++;
}
}
let H1 = new Array();
for (let i1 = 0; i1 < n; i1++)
H1[i1] = new Array();
let H2 = new Array();
for (let i1 = 0; i1 < n; i1++)
H2[i1] = new Array();
ind = frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2);
// 出力
if (ind > 0)
document.getElementById("ans").value = "解を求めることができません！";
else {
let str = "";
for (let i1 = 0; i1 < n; i1++)
str = str + rl[i1] + " i " + im[i1] + "\n";
document.getElementById("ans").value = str;
}
}
/*************************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b, c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* A : 行列 */
/* H1, H2 : 作業域（nxnの行列） */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
function frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2)
{
let i1;
let i2;
let i3;
let i4;
let ind;
a[0] = 1.0;
// b1の計算
a[1] = 0.0;
for (i1 = 0; i1 < n; i1++)
a[1] -= A[i1][i1];
// b2の計算
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
H1[i1][i2] = A[i1][i2];
H1[i1][i1] += a[1];
}
a[2] = 0.0;
for (i1 = 0; i1 < n; i1++) {
for (i2 = 0; i2 < n; i2++)
a[2] -= A[i1][i2] * H1[i2][i1];
}
a[2] *= 0.5;
// b3からbnの計算
for (i1 = 3; i1 <= n; i1++) {
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++) {
H2[i2][i3] = 0.0;
for (i4 = 0; i4 < n; i4++)
H2[i2][i3] += A[i2][i4] * H1[i4][i3];
}
H2[i2][i2] += a[i1-1];
}
a[i1] = 0.0;
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++)
a[i1] -= A[i2][i3] * H2[i3][i2];
}
a[i1] /= i1;
for (i2 = 0; i2 < n; i2++) {
for (i3 = 0; i3 < n; i3++)
H1[i2][i3] = H2[i2][i3];
}
}
// ベアストウ法
ind = bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);
return ind;
}
/*************************************************/
/* 実係数代数方程式の解（ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b,c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* p : 答えの位置 */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
function bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p)
{
let D;
let dp;
let dq;
let p1 = p0;
let p2 = 0.0;
let q1 = q0;
let q2 = 0.0;
let i1;
let ind = 0;
let count = 0;
/*
１次の場合
*/
if (n == 1) {
if (Math.abs(a[0]) < eps)
ind = 1;
else {
rl[p] = -a[1] / a[0];
im[p] = 0.0;
}
}
/*
２次の場合
*/
else if (n == 2) {
// １次式
if (Math.abs(a[0]) < eps) {
if (Math.abs(a[1]) < eps)
ind = 1;
else {
rl[p] = -a[2] / a[1];
im[p] = 0.0;
}
}
// ２次式
else {
D = a[1] * a[1] - 4.0 * a[0] * a[2];
if (D < 0.0) { // 虚数
D = Math.sqrt(-D);
a[0] *= 2.0;
rl[p] = -a[1] / a[0];
rl[p+1] = -a[1] / a[0];
im[p] = D / a[0];
im[p+1] = -im[p];
}
else { // 実数
D = Math.sqrt(D);
a[0] = 1.0 / (2.0 * a[0]);
rl[p] = a[0] * (-a[1] + D);
rl[p+1] = a[0] * (-a[1] - D);
im[p] = 0.0;
im[p+1] = 0.0;
}
}
}
// ３次以上の場合
else {
// 因数分解
ind = 1;
while (ind > 0 && count <= ct) {
for (i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
b[i1] = a[i1];
else if (i1 == 1)
b[i1] = a[i1] - p1 * b[i1-1];
else
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2];
}
for (i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
c[i1] = b[i1];
else if (i1 == 1)
c[i1] = b[i1] - p1 * c[i1-1];
else
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2];
}
D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1]);
if (Math.abs(D) < eps)
return ind;
else {
dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / D;
dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / D;
p2 = p1 + dp;
q2 = q1 + dq;
if (Math.abs(dp) < eps && Math.abs(dq) < eps)
ind = 0;
else {
count++;
p1 = p2;
q1 = q2;
}
}
}
if (ind == 0) {
// ２次方程式を解く
D = p2 * p2 - 4.0 * q2;
if (D < 0.0) { // 虚数
D = Math.sqrt(-D);
rl[p] = -0.5 * p2;
rl[p+1] = -0.5 * p2;
im[p] = 0.5 * D;
im[p+1] = -im[p];
}
else { // 実数
D = Math.sqrt(D);
rl[p] = 0.5 * (-p2 + D);
rl[p+1] = 0.5 * (-p2 - D);
im[p] = 0.0;
im[p+1] = 0.0;
}
// 残りの方程式を解く
n -= 2;
for (i1 = 0; i1 <= n; i1++)
a[i1] = b[i1];
ind = bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p+2);
}
}
return ind;
}
</SCRIPT>
</HEAD>
<BODY STYLE="font-size: 130%; background-color: #eeffee;">
<H2 STYLE="text-align:center"><B>行列の固有値（フレーム法＋ベアストウ法）</B></H2>
<DL>
<DT> テキストフィールドおよびテキストエリアには，例として，以下に示す行列の固有値を求める場合に対する値が設定されています．他の問題を実行する場合は，それらを適切に修正してください．
<P STYLE="text-align:center"><IMG SRC="eigen.gif"></P>
</DL>
<DIV STYLE="text-align:center">
次数：<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="3">
p0：<INPUT ID="p0" STYLE="font-size: 100%;" TYPE="text" SIZE="2" VALUE="0">
q0：<INPUT ID="q0" STYLE="font-size: 100%;" TYPE="text" SIZE="2" VALUE="0">
最大繰り返し回数：<INPUT ID="trial" STYLE="font-size: 100%;" TYPE="text" SIZE="4" VALUE="1000">
<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
<TEXTAREA ID="ar" COLS="50" ROWS="15" STYLE="font-size: 100%">7 2 1
2 1 -4
1 -4 2</TEXTAREA><BR><BR>
<TEXTAREA ID="ans" COLS="50" ROWS="15" STYLE="font-size: 100%"></TEXTAREA>
</DIV>
</BODY>
</HTML>
<?php
/********************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* coded by Y.Suganuma */
/********************************************/
// データの設定
$ct = 1000;
$eps = 1.0e-10;
$p0 = 0.0;
$q0 = 0.0;
$n = 3;
$a = array($n+1);
$b = array($n+1);
$c = array($n+1);
$rl = array($n);
$im = array($n);
$A = array($n);
$H1 = array($n);
$H2 = array($n);
for ($i1 = 0; $i1 < $n; $i1++) {
$A[$i1] = array($n);
$H1[$i1] = array($n);
$H2[$i1] = array($n);
}
$A[0][0] = 7.0;
$A[0][1] = 2.0;
$A[0][2] = 1.0;
$A[1][0] = 2.0;
$A[1][1] = 1.0;
$A[1][2] = -4.0;
$A[2][0] = 1.0;
$A[2][1] = -4.0;
$A[2][2] = 2.0;
// 計算
$ind = Frame($n, $ct, $eps, $p0, $q0, $a, $b, $c, $rl, $im, $A, $H1, $H2);
// 出力
if ($ind > 0)
printf("収束しませんでした！\n");
else {
for ($i1 = 0; $i1 < $n; $i1++)
printf(" %f i %f\n", $rl[$i1], $im[$i1]);
}
/*************************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b, c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* A : 行列 */
/* H1, H2 : 作業域（nxnの行列） */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
function Frame($n, $ct, $eps, $p0, $q0, $a, $b, $c, &$rl, &$im, $A, $H1, $H2)
{
$a[0] = 1.0;
// a1の計算
$a[1] = 0.0;
for ($i1 = 0; $i1 < $n; $i1++)
$a[1] -= $A[$i1][$i1];
// a2の計算
for ($i1 = 0; $i1 < $n; $i1++) {
for ($i2 = 0; $i2 < $n; $i2++)
$H1[$i1][$i2] = $A[$i1][$i2];
$H1[$i1][$i1] += $a[1];
}
$a[2] = 0.0;
for ($i1 = 0; $i1 < $n; $i1++) {
for ($i2 = 0; $i2 < $n; $i2++)
$a[2] -= $A[$i1][$i2] * $H1[$i2][$i1];
}
$a[2] *= 0.5;
// a3からanの計算
for ($i1 = 3; $i1 <= $n; $i1++) {
for ($i2 = 0; $i2 < $n; $i2++) {
for ($i3 = 0; $i3 < $n; $i3++) {
$H2[$i2][$i3] = 0.0;
for ($i4 = 0; $i4 < $n; $i4++)
$H2[$i2][$i3] += $A[$i2][$i4] * $H1[$i4][$i3];
}
$H2[$i2][$i2] += $a[$i1-1];
}
$a[$i1] = 0.0;
for ($i2 = 0; $i2 < $n; $i2++) {
for ($i3 = 0; $i3 < $n; $i3++)
$a[$i1] -= $A[$i2][$i3] * $H2[$i3][$i2];
}
$a[$i1] /= $i1;
for ($i2 = 0; $i2 < $n; $i2++) {
for ($i3 = 0; $i3 < $n; $i3++)
$H1[$i2][$i3] = $H2[$i2][$i3];
}
}
// ベアストウ法
$ind = Bairstow($n, $ct, $eps, $p0, $q0, $a, $b, $c, $rl, $im, 0);
return $ind;
}
/*************************************************/
/* 実係数代数方程式の解（ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b,c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* k : 結果を設定する配列の位置 */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
function Bairstow($n, $ct, $eps, $p0, $q0, $a, $b, $c, &$rl, &$im, $k)
{
$p1 = $p0;
$p2 = 0.0;
$q1 = $q0;
$q2 = 0.0;
$ind = 0;
$count = 0;
/*
１次の場合
*/
if ($n == 1) {
if (abs($a[0]) < $eps)
$ind = 1;
else {
$rl[$k] = -$a[1] / $a[0];
$im[$k] = 0.0;
}
}
/*
２次の場合
*/
else if ($n == 2) {
// １次式
if (abs($a[0]) < $eps) {
if (abs($a[1]) < $eps)
$ind = 1;
else {
$rl[$k] = -$a[2] / $a[1];
$im[$k] = 0.0;
}
}
// ２次式
else {
$D = $a[1] * $a[1] - 4.0 * $a[0] * $a[2];
if ($D < 0.0) { // 虚数
$D = sqrt(-$D);
$a[0] *= 2.0;
$rl[$k] = -$a[1] / $a[0];
$rl[$k+1] = -$a[1] / $a[0];
$im[$k] = $D / $a[0];
$im[$k+1] = -$im[$k];
}
else { // 実数
$D = sqrt($D);
$a[0] = 1.0 / (2.0 * $a[0]);
$rl[$k] = $a[0] * (-$a[1] + $D);
$rl[$k+1] = $a[0] * (-$a[1] - $D);
$im[$k] = 0.0;
$im[$k+1] = 0.0;
}
}
}
// ３次以上の場合
else {
// 因数分解
$ind = 1;
while ($ind > 0 && $count <= $ct) {
for ($i1 = 0; $i1 <= $n; $i1++) {
if ($i1 == 0)
$b[$i1] = $a[$i1];
else if ($i1 == 1)
$b[$i1] = $a[$i1] - $p1 * $b[$i1-1];
else
$b[$i1] = $a[$i1] - $p1 * $b[$i1-1] - $q1 * $b[$i1-2];
}
for ($i1 = 0; $i1 <= $n; $i1++) {
if ($i1 == 0)
$c[$i1] = $b[$i1];
else if ($i1 == 1)
$c[$i1] = $b[$i1] - $p1 * $c[$i1-1];
else
$c[$i1] = $b[$i1] - $p1 * $c[$i1-1] - $q1 * $c[$i1-2];
}
$D = $c[$n-2] * $c[$n-2] - $c[$n-3] * ($c[$n-1] - $b[$n-1]);
if (abs($D) < $eps)
return $ind;
else {
$dp = ($b[$n-1] * $c[$n-2] - $b[$n] * $c[$n-3]) / $D;
$dq = ($b[$n] * $c[$n-2] - $b[$n-1] * ($c[$n-1] - $b[$n-1])) / $D;
$p2 = $p1 + $dp;
$q2 = $q1 + $dq;
if (abs($dp) < $eps && abs($dq) < $eps)
$ind = 0;
else {
$count++;
$p1 = $p2;
$q1 = $q2;
}
}
}
if ($ind == 0) {
// ２次方程式を解く
$D = $p2 * $p2 - 4.0 * $q2;
if ($D < 0.0) { // 虚数
$D = sqrt(-$D);
$rl[$k] = -0.5 * $p2;
$rl[$k+1] = -0.5 * $p2;
$im[$k] = 0.5 * $D;
$im[$k+1] = -$im[$k];
}
else { // 実数
$D = sqrt($D);
$rl[$k] = 0.5 * (-$p2 + $D);
$rl[$k+1] = 0.5 * (-$p2 - $D);
$im[$k] = 0.0;
$im[$k+1] = 0.0;
}
// 残りの方程式を解く
$n -= 2;
for ($i1 = 0; $i1 <= $n; $i1++)
$a[$i1] = $b[$i1];
$ind = Bairstow($n, $ct, $eps, $p0, $q0, $a, $b, $c, $rl, $im, $k+2);
}
}
return $ind;
}
?>
#*******************************************/
# 行列の固有値（フレーム法＋ベアストウ法） */
# coded by Y.Suganuma */
#*******************************************/
#************************************************/
# 行列の固有値（フレーム法＋ベアストウ法） */
# n : 次数 */
# ct : 最大繰り返し回数 */
# eps : 収束判定条件 */
# p0, q0 : x2+px+qにおけるp,qの初期値 */
# a : 係数（最高次から与え，値は変化する） */
# b, c : 作業域（(n+1)次の配列） */
# rl, im : 結果の実部と虚部 */
# aa : 行列 */
# h1, h2 : 作業域（nxnの行列） */
# return : =0 : 正常 */
# =1 : 収束せず */
# coded by Y.Suganuma */
#************************************************/
def Frame(n, ct, eps, p0, q0, a, b, c, rl, im, aa, h1, h2)
a[0] = 1.0
# a1の計算
a[1] = 0.0
for i1 in 0 ... n
a[1] -= aa[i1][i1]
end
# a2の計算
for i1 in 0 ... n
for i2 in 0 ... n
h1[i1][i2] = aa[i1][i2]
end
h1[i1][i1] += a[1]
end
a[2] = 0.0
for i1 in 0 ... n
for i2 in 0 ... n
a[2] -= aa[i1][i2] * h1[i2][i1]
end
end
a[2] *= 0.5
# a3からanの計算
for i1 in 3 ... n+1
for i2 in 0 ... n
for i3 in 0 ... n
h2[i2][i3] = 0.0
for i4 in 0 ... n
h2[i2][i3] += aa[i2][i4] * h1[i4][i3]
end
end
h2[i2][i2] += a[i1-1]
end
a[i1] = 0.0
for i2 in 0 ... n
for i3 in 0 ... n
a[i1] -= aa[i2][i3] * h2[i3][i2]
end
end
a[i1] /= i1
for i2 in 0 ... n
for i3 in 0 ... n
h1[i2][i3] = h2[i2][i3]
end
end
end
# ベアストウ法
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
return ind
end
#************************************************/
# 実係数代数方程式の解（ベアストウ法） */
# n : 次数 */
# ct : 最大繰り返し回数 */
# eps : 収束判定条件 */
# p0, q0 : x2+px+qにおけるp,qの初期値 */
# a : 係数（最高次から与え，値は変化する） */
# b,c : 作業域（(n+1)次の配列） */
# rl, im : 結果の実部と虚部 */
# k : 結果の位置 */
# return : =0 : 正常 */
# =1 : 収束せず */
# coded by Y.Suganuma */
#************************************************/
def Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k)
# 初期設定
p1 = p0
p2 = 0.0
q1 = q0
q2 = 0.0
ind = 0
count = 0
#
# １次の場合
#
if n == 1
if a[0].abs() < eps
ind = 1
else
rl[k] = -a[1] / a[0]
im[k] = 0.0
end
#
# ２次の場合
#
elsif n == 2
# １次式
if a[0].abs() < eps
if a[1].abs() < eps
ind = 1
else
rl[k] = -a[2] / a[1]
im[k] = 0.0
end
# ２次式
else
d = a[1] * a[1] - 4.0 * a[0] * a[2]
if d < 0.0 # 虚数
d = Math.sqrt(-d)
a[0] *= 2.0
rl[k] = -a[1] / a[0]
rl[k+1] = -a[1] / a[0]
im[k] = d / a[0]
im[k+1] = -im[0]
else # 実数
d = Math.sqrt(d)
a[0] = 1.0 / (2.0 * a[0])
rl[k] = a[0] * (-a[1] + d)
rl[k+1] = a[0] * (-a[1] - d)
im[k] = 0.0
im[k+1] = 0.0
end
end
# ３次以上の場合
else
# 因数分解
ind = 1
while ind > 0 && count <= ct
for i1 in 0 ... n+1
if i1 == 0
b[i1] = a[i1]
elsif i1 == 1
b[i1] = a[i1] - p1 * b[i1-1]
else
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2]
end
end
for i1 in 0 ... n+1
if i1 == 0
c[i1] = b[i1]
elsif i1 == 1
c[i1] = b[i1] - p1 * c[i1-1]
else
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2]
end
end
d = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1])
if d.abs() < eps
return ind
else
dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / d
dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / d
p2 = p1 + dp
q2 = q1 + dq
if dp.abs() < eps && dq.abs() < eps
ind = 0
else
count += 1
p1 = p2
q1 = q2
end
end
end
if ind == 0
# ２次方程式を解く
d = p2 * p2 - 4.0 * q2
if d < 0.0 # 虚数
d = Math.sqrt(-d)
rl[k] = -0.5 * p2
rl[k+1] = -0.5 * p2
im[k] = 0.5 * d
im[k+1] = -im[k]
else # 実数
d = Math.sqrt(d)
rl[k] = 0.5 * (-p2 + d)
rl[k+1] = 0.5 * (-p2 - d)
im[k] = 0.0
im[k+1] = 0.0
end
# 残りの方程式を解く
n -= 2
for i1 in 0 ... n+1
a[i1] = b[i1]
end
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k+2)
end
end
return ind
end
# データの設定
ct = 1000
eps = 1.0e-10
p0 = 0.0
q0 = 0.0
n = 3
a = Array.new(n+1)
b = Array.new(n+1)
c = Array.new(n+1)
rl = Array.new(n)
im = Array.new(n)
aa = Array.new(n)
h1 = Array.new(n)
h2 = Array.new(n)
for i1 in 0 ... n
aa[i1] = Array.new(n)
h1[i1] = Array.new(n)
h2[i1] = Array.new(n)
end
aa[0][0] = 7.0
aa[0][1] = 2.0
aa[0][2] = 1.0
aa[1][0] = 2.0
aa[1][1] = 1.0
aa[1][2] = -4.0
aa[2][0] = 1.0
aa[2][1] = -4.0
aa[2][2] = 2.0
# 計算
ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, aa, h1, h2)
# 出力
if ind > 0
printf("収束しませんでした！\n")
else
for i1 in 0 ... n
printf(" %f i %f\n", rl[i1], im[i1])
end
end
# -*- coding: UTF-8 -*-
import numpy as np
from math import *
############################################
# 行列の固有値（フレーム法＋ベアストウ法）
# n : 次数
# ct : 最大繰り返し回数
# eps : 収束判定条件
# p0, q0 : x2+px+qにおけるp,qの初期値
# a : 係数（最高次から与え，値は変化する）
# b, c : 作業域（(n+1)次の配列）
# r : 結果
# A : 行列
# H1, H2 : 作業域（nxnの行列）
# return : =0 : 正常
# =1 : 収束せず
# coded by Y.Suganuma
############################################
def Frame(n, ct, eps, p0, q0, a, b, c, r, A, H1, H2) :
a[0] = 1.0
# a1の計算
a[1] = 0.0
for i1 in range(0, n) :
a[1] -= A[i1][i1]
# a2の計算
for i1 in range(0, n) :
for i2 in range(0, n) :
H1[i1][i2] = A[i1][i2]
H1[i1][i1] += a[1]
a[2] = 0.0
for i1 in range(0, n) :
for i2 in range(0, n) :
a[2] -= A[i1][i2] * H1[i2][i1]
a[2] *= 0.5
# a3からanの計算
for i1 in range(3, n+1) :
for i2 in range(0, n) :
for i3 in range(0, n) :
H2[i2][i3] = 0.0
for i4 in range(0, n) :
H2[i2][i3] += A[i2][i4] * H1[i4][i3]
H2[i2][i2] += a[i1-1]
a[i1] = 0.0
for i2 in range(0, n) :
for i3 in range(0, n) :
a[i1] -= A[i2][i3] * H2[i3][i2]
a[i1] /= i1
for i2 in range(0, n) :
for i3 in range(0, n) :
H1[i2][i3] = H2[i2][i3]
# ベアストウ法
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, 0)
return ind
############################################
# 実係数代数方程式の解（ベアストウ法）
# n : 次数
# ct : 最大繰り返し回数
# eps : 収束判定条件
# p0, q0 : x2+px+qにおけるp,qの初期値
# a : 係数（最高次から与え，値は変化する）
# b,c : 作業域（(n+1)次の配列）
# r : 結果
# k : 結果の位置
# return : =0 : 正常
# =1 : 収束せず
# coded by Y.Suganuma
############################################
def Bairstow(n, ct, eps, p0, q0, a, b, c, r, k) :
p1 = p0
p2 = 0.0
q1 = q0
q2 = 0.0
ind = 0
count = 0
# １次の場合
if n == 1 :
if abs(a[0]) < eps :
ind = 1
else :
r[k] = complex(-a[1] / a[0], 0)
# ２次の場合
elif n == 2 :
# １次式
if abs(a[0]) < eps :
if abs(a[1]) < eps :
ind = 1
else :
r[k] = complex(-a[2] / a[1], 0)
# ２次式
else :
D = a[1] * a[1] - 4.0 * a[0] * a[2]
if D < 0.0 : # 虚数
D = sqrt(-D)
a[0] *= 2.0
r[k] = complex(-a[1] / a[0], D / a[0])
r[k+1] = complex(-a[1] / a[0], -D / a[0])
else : # 実数
D = sqrt(D)
a[0] = 1.0 / (2.0 * a[0])
r[k] = complex(a[0] * (-a[1] + D), 0)
r[k+1] = complex(a[0] * (-a[1] - D), 0)
# ３次以上の場合
else :
# 因数分解
ind = 1
while ind > 0 and count <= ct :
for i1 in range(0, n+1) :
if i1 == 0 :
b[i1] = a[i1]
elif i1 == 1 :
b[i1] = a[i1] - p1 * b[i1-1]
else :
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2]
for i1 in range(0, n+1) :
if i1 == 0 :
c[i1] = b[i1]
elif i1 == 1 :
c[i1] = b[i1] - p1 * c[i1-1]
else :
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2]
D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1])
if fabs(D) < eps :
return ind
else :
dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / D
dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / D
p2 = p1 + dp
q2 = q1 + dq
if abs(dp) < eps and fabs(dq) < eps :
ind = 0
else :
count += 1
p1 = p2
q1 = q2
if ind == 0 :
# ２次方程式を解く
D = p2 * p2 - 4.0 * q2
if D < 0.0 : # 虚数
D = sqrt(-D)
r[k] = complex(-0.5 * p2, 0.5 * D)
r[k+1] = complex(-0.5 * p2, -0.5 * D)
else : # 実数
D = sqrt(D)
r[k] = complex(0.5 * (-p2 + D), 0)
r[k+1] = complex(0.5 * (-p2 - D), 0)
# 残りの方程式を解く
n -= 2
for i1 in range(0, n+1) :
a[i1] = b[i1]
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, k+2)
return ind
############################################
# 行列の固有値（フレーム法＋ベアストウ法） */
# coded by Y.Suganuma */
############################################
# データの設定
ct = 1000
eps = 1.0e-10
p0 = 0.0
q0 = 0.0
n = 3
a = np.empty(n+1, np.float)
b = np.empty(n+1, np.float)
c = np.empty(n+1, np.float)
r = np.empty(n, np.complex)
A = np.array([[7, 2, 1],[2, 1, -4],[1, -4, 2]], np.float)
H1 = np.empty((n, n), np.float)
H2 = np.empty((n, n), np.float)
# 計算
ind = Frame(n, ct, eps, p0, q0, a, b, c, r, A, H1, H2)
# 出力
if ind > 0 :
print("収束しませんでした！")
else :
for i1 in range(0, n) :
print(" " + str(r[i1]))
/********************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* coded by Y.Suganuma */
/********************************************/
using System;
class Program
{
static void Main()
{
Test1 ts = new Test1();
}
}
class Test1
{
public Test1()
{
// データの設定
int ct = 1000;
int n = 3;
double eps = 1.0e-10;
double p0 = 0.0;
double q0 = 0.0;
double[] a = new double [n+1];
double[] b = new double [n+1];
double[] c = new double [n+1];
double[] rl = new double [n];
double[] im = new double [n];
double[][] A = new double [][] {
new double[] {7.0, 2.0, 1.0},
new double[] {2.0, 1.0, -4.0},
new double[] {1.0, -4.0, 2.0}
};
double[][] H1 = new double [n][];
double[][] H2 = new double [n][];
for (int i1 = 0; i1 < n; i1++) {
H1[i1] = new double [n];
H2[i1] = new double [n];
}
// 計算
int ind = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, A, H1, H2);
// 出力
if (ind > 0)
Console.WriteLine("収束しませんでした！");
else {
for (int i1 = 0; i1 < n; i1++)
Console.WriteLine(" " + rl[i1] + " i " + im[i1]);
}
}
/*************************************************/
/* 行列の固有値（フレーム法＋ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b, c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* A : 行列 */
/* H1, H2 : 作業域（nxnの行列） */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
int Frame(int n, int ct, double eps, double p0, double q0, double[] a, double[] b,
double[] c, double[] rl, double[] im, double[][] A, double[][] H1,
double[][] H2)
{
a[0] = 1.0;
// a1の計算
a[1] = 0.0;
for (int i1 = 0; i1 < n; i1++)
a[1] -= A[i1][i1];
// a2の計算
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
H1[i1][i2] = A[i1][i2];
H1[i1][i1] += a[1];
}
a[2] = 0.0;
for (int i1 = 0; i1 < n; i1++) {
for (int i2 = 0; i2 < n; i2++)
a[2] -= A[i1][i2] * H1[i2][i1];
}
a[2] *= 0.5;
// a3からanの計算
for (int i1 = 3; i1 <= n; i1++) {
for (int i2 = 0; i2 < n; i2++) {
for (int i3 = 0; i3 < n; i3++) {
H2[i2][i3] = 0.0;
for (int i4 = 0; i4 < n; i4++)
H2[i2][i3] += A[i2][i4] * H1[i4][i3];
}
H2[i2][i2] += a[i1-1];
}
a[i1] = 0.0;
for (int i2 = 0; i2 < n; i2++) {
for (int i3 = 0; i3 < n; i3++)
a[i1] -= A[i2][i3] * H2[i3][i2];
}
a[i1] /= i1;
for (int i2 = 0; i2 < n; i2++) {
for (int i3 = 0; i3 < n; i3++)
H1[i2][i3] = H2[i2][i3];
}
}
// ベアストウ法
int ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0);
return ind;
}
/*************************************************/
/* 実係数代数方程式の解（ベアストウ法） */
/* n : 次数 */
/* ct : 最大繰り返し回数 */
/* eps : 収束判定条件 */
/* p0, q0 : x2+px+qにおけるp,qの初期値 */
/* a : 係数（最高次から与え，値は変化する） */
/* b,c : 作業域（(n+1)次の配列） */
/* rl, im : 結果の実部と虚部 */
/* p : 答えの位置 */
/* return : =0 : 正常 */
/* =1 : 収束せず */
/* coded by Y.Suganuma */
/*************************************************/
static int Bairstow(int n, int ct, double eps, double p0, double q0, double[] a,
double[] b, double[] c, double[] rl, double[] im, int p)
{
int ind = 0;
/*
１次の場合
*/
if (n == 1) {
if (Math.Abs(a[0]) < eps)
ind = 1;
else {
rl[p] = -a[1] / a[0];
im[p] = 0.0;
}
}
/*
２次の場合
*/
else if (n == 2) {
double D;
// １次式
if (Math.Abs(a[0]) < eps) {
if (Math.Abs(a[1]) < eps)
ind = 1;
else {
rl[p] = -a[2] / a[1];
im[p] = 0.0;
}
}
// ２次式
else {
D = a[1] * a[1] - 4.0 * a[0] * a[2];
if (D < 0.0) { // 虚数
D = Math.Sqrt(-D);
a[0] *= 2.0;
rl[p] = -a[1] / a[0];
rl[p+1] = -a[1] / a[0];
im[p] = D / a[0];
im[p+1] = -im[p];
}
else { // 実数
D = Math.Sqrt(D);
a[0] = 1.0 / (2.0 * a[0]);
rl[p] = a[0] * (-a[1] + D);
rl[p+1] = a[0] * (-a[1] - D);
im[p] = 0.0;
im[p+1] = 0.0;
}
}
}
// ３次以上の場合
else {
// 因数分解
ind = 1;
int count = 0;
double D, dp, dq, p1 = p0, p2 = 0.0, q1 = q0, q2 = 0.0;
while (ind > 0 && count <= ct) {
for (int i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
b[i1] = a[i1];
else if (i1 == 1)
b[i1] = a[i1] - p1 * b[i1-1];
else
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2];
}
for (int i1 = 0; i1 <= n; i1++) {
if (i1 == 0)
c[i1] = b[i1];
else if (i1 == 1)
c[i1] = b[i1] - p1 * c[i1-1];
else
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2];
}
D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1]);
if (Math.Abs(D) < eps)
return ind;
else {
dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / D;
dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / D;
p2 = p1 + dp;
q2 = q1 + dq;
if (Math.Abs(dp) < eps && Math.Abs(dq) < eps)
ind = 0;
else {
count++;
p1 = p2;
q1 = q2;
}
}
}
if (ind == 0) {
// ２次方程式を解く
D = p2 * p2 - 4.0 * q2;
if (D < 0.0) { // 虚数
D = Math.Sqrt(-D);
rl[p] = -0.5 * p2;
rl[p+1] = -0.5 * p2;
im[p] = 0.5 * D;
im[p+1] = -im[p];
}
else { // 実数
D = Math.Sqrt(D);
rl[p] = 0.5 * (-p2 + D);
rl[p+1] = 0.5 * (-p2 - D);
im[p] = 0.0;
im[p+1] = 0.0;
}
// 残りの方程式を解く
n -= 2;
for (int i1 = 0; i1 <= n; i1++)
a[i1] = b[i1];
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p+2);
}
}
return ind;
}
}
''''''''''''''''''''''''''''''''''''''''''''
' 行列の固有値（フレーム法＋ベアストウ法） '
' coded by Y.Suganuma '
''''''''''''''''''''''''''''''''''''''''''''
Module Test
Sub Main()
' データの設定
Dim ct As Integer = 1000
Dim n As Integer = 3
Dim eps As Double = 1.0e-10
Dim p0 As Double = 0.0
Dim q0 As Double = 0.0
Dim a(n+1) As Double
Dim b(n+1) As Double
Dim c(n+1) As Double
Dim rl(n) As Double
Dim im(n) As Double
Dim AA(,) As Double = {{7.0, 2.0, 1.0}, {2.0, 1.0, -4.0}, {1.0, -4.0, 2.0}}
Dim H1(n,n) As Double
Dim H2(n,n) As Double
' 計算
Dim ind As Integer = Frame(n, ct, eps, p0, q0, a, b, c, rl, im, AA, H1, H2)
' 出力
If ind > 0
Console.WriteLine("収束しませんでした！")
Else
For i1 As Integer = 0 To n-1
Console.WriteLine(" " & rl(i1) & " i " & im(i1))
Next
End If
End Sub
'''''''''''''''''''''''''''''''''''''''''''''''''
' 行列の固有値（フレーム法＋ベアストウ法） '
' n : 次数 '
' ct : 最大繰り返し回数 '
' eps : 収束判定条件 '
' p0, q0 : x2+px+qにおけるp,qの初期値 '
' a : 係数（最高次から与え，値は変化する） '
' b, c : 作業域（(n+1)次の配列） '
' rl, im : 結果の実部と虚部 '
' AA : 行列 '
' H1, H2 : 作業域（nxnの行列） '
' return : =0 : 正常 '
' =1 : 収束せず '
' coded by Y.Suganuma '
'''''''''''''''''''''''''''''''''''''''''''''''''
Function Frame(n As Integer, ct As Integer, eps As Double, p0 As Double,
q0 As Double, a() As Double, b() As Double, c() As Double,
rl() As Double, im() As Double, AA(,) As Double, H1(,) As Double,
H2(,) As Double)
a(0) = 1.0
' a1の計算
a(1) = 0.0
For i1 As Integer = 0 To n-1
a(1) -= AA(i1,i1)
Next
' a2の計算
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
H1(i1,i2) = AA(i1,i2)
Next
H1(i1,i1) += a(1)
Next
a(2) = 0.0
For i1 As Integer = 0 To n-1
For i2 As Integer = 0 To n-1
a(2) -= AA(i1,i2) * H1(i2,i1)
Next
Next
a(2) *= 0.5
' a3からanの計算
For i1 As Integer = 3 To n
For i2 As Integer = 0 To n-1
For i3 As Integer = 0 To n-1
H2(i2,i3) = 0.0
For i4 As Integer = 0 To n-1
H2(i2,i3) += AA(i2,i4) * H1(i4,i3)
Next
Next
H2(i2,i2) += a(i1-1)
Next
a(i1) = 0.0
For i2 As Integer = 0 To n-1
For i3 As Integer = 0 To n-1
a(i1) -= AA(i2,i3) * H2(i3,i2)
Next
Next
a(i1) /= i1
For i2 As Integer = 0 To n-1
For i3 As Integer = 0 To n-1
H1(i2,i3) = H2(i2,i3)
Next
Next
Next
' ベアストウ法
Dim ind As Integer = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
Return ind
End Function
'''''''''''''''''''''''''''''''''''''''''''''''''
' 実係数代数方程式の解（ベアストウ法） '
' n : 次数 '
' ct : 最大繰り返し回数 '
' eps : 収束判定条件 '
' p0, q0 : x2+px+qにおけるp,qの初期値 '
' a : 係数（最高次から与え，値は変化する） '
' b,c : 作業域（(n+1)次の配列） '
' rl, im : 結果の実部と虚部 '
' p : 答えの位置 '
' return : =0 : 正常 '
' =1 : 収束せず '
' coded by Y.Suganuma '
'''''''''''''''''''''''''''''''''''''''''''''''''
Function Bairstow(n As Integer, ct As Integer, eps As Double, p0 As Double,
q0 As Double, a() As Double, b() As Double, c() As Double,
rl() As Double, im() As Double, p As Integer)
Dim ind As Integer = 0
'
' １次の場合
'
If n = 1
If Math.Abs(a(0)) < eps
ind = 1
Else
rl(p) = -a(1) / a(0)
im(p) = 0.0
End If
'
' ２次の場合
'
ElseIf n = 2
Dim D As Double
' １次式
If Math.Abs(a(0)) < eps
If Math.Abs(a(1)) < eps
ind = 1
Else
rl(p) = -a(2) / a(1)
im(p) = 0.0
End If
' ２次式
Else
D = a(1) * a(1) - 4.0 * a(0) * a(2)
If D < 0.0 ' 虚数
D = Math.Sqrt(-D)
a(0) *= 2.0
rl(p) = -a(1) / a(0)
rl(p+1) = -a(1) / a(0)
im(p) = D / a(0)
im(p+1) = -im(p)
Else ' 実数
D = Math.Sqrt(D)
a(0) = 1.0 / (2.0 * a(0))
rl(p) = a(0) * (-a(1) + D)
rl(p+1) = a(0) * (-a(1) - D)
im(p) = 0.0
im(p+1) = 0.0
End If
End If
' ３次以上の場合
Else
' 因数分解
ind = 1
Dim count As Integer = 0
Dim D As Double
Dim dp As Double
Dim dq As Double
Dim p1 As Double = p0
Dim p2 As Double = 0.0
Dim q1 As Double = q0
Dim q2 As Double = 0.0
Do While ind > 0 and count <= ct
For i1 As Integer = 0 To n
If i1 = 0
b(i1) = a(i1)
ElseIf i1 = 1
b(i1) = a(i1) - p1 * b(i1-1)
Else
b(i1) = a(i1) - p1 * b(i1-1) - q1 * b(i1-2)
End If
Next
For i1 As Integer = 0 To n
If i1 = 0
c(i1) = b(i1)
ElseIf i1 = 1
c(i1) = b(i1) - p1 * c(i1-1)
Else
c(i1) = b(i1) - p1 * c(i1-1) - q1 * c(i1-2)
End If
Next
D = c(n-2) * c(n-2) - c(n-3) * (c(n-1) - b(n-1))
If Math.Abs(D) < eps
Return ind
Else
dp = (b(n-1) * c(n-2) - b(n) * c(n-3)) / D
dq = (b(n) * c(n-2) - b(n-1) * (c(n-1) - b(n-1))) / D
p2 = p1 + dp
q2 = q1 + dq
If Math.Abs(dp) < eps and Math.Abs(dq) < eps
ind = 0
Else
count += 1
p1 = p2
q1 = q2
End If
End If
Loop
If ind = 0
' ２次方程式を解く
D = p2 * p2 - 4.0 * q2
If D < 0.0 ' 虚数
D = Math.Sqrt(-D)
rl(p) = -0.5 * p2
rl(p+1) = -0.5 * p2
im(p) = 0.5 * D
im(p+1) = -im(p)
Else ' 実数
D = Math.Sqrt(D)
rl(p) = 0.5 * (-p2 + D)
rl(p+1) = 0.5 * (-p2 - D)
im(p) = 0.0
im(p+1) = 0.0
End If
' 残りの方程式を解く
n -= 2
For i1 As Integer = 0 To n
a(i1) = b(i1)
Next
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, p+2)
End If
End If
Return ind
End Function
End Module
| 情報学部 | 菅沼ホーム | 目次 | 索引 |