|
投稿者：tibita
投稿日：2014年10月14日(火)17時32分32秒
|
|
|
|
C言語のプログラミングについて勉強しています。
8次の代数方程式をDK法で解き、実部、虚部の解を共にモニタ上に出力するプログラムを作りたいのですが、C言語で書かれたサンプルプログラムを書いていただけませんか？
|
|
|
|
投稿者：山中和義
投稿日：2014年10月16日(木)10時37分3秒
|
|
|
|
> No.3525[元記事へ]
tibitaさんへのお返事です。
> C言語のプログラミングについて勉強しています。
> 8次の代数方程式をDK法で解き、実部、虚部の解を共にモニタ上に出力するプログラムを作りたいのですが、C言語で書かれたサンプルプログラムを書いていただけませんか？
代数方程式
p(z)=z^8-36z^7+546z^6-4536z^5+22449z^4-67284z^3+118124z^2-109584z+40320
の解
最近の言語では、COMPLEX（複素数）を扱えるので、記述はより簡潔になると思います。
下記のANSI C程度のサンプルは、実部と虚部（ x1(n), x2(n) ）を分離して計算しています。
//DKA（Durand-Kerner-Aberth）法によるn次代数方程式f(x)=0の解法
#include <stdio.h>
#include <math.h>
#define PI 3.14159265358979323846
#define N 8 //n次
int main(void)
{
double a[N+1] = {40320,-109584,118124,-67284,22449,-4536,546,-36,1}; //係数
double x1[N],x2[N]; //n個の根（実部、虚部）
int i,j;
for (i=0; i<=N; i++) a[i]/=a[N]; //x^n+a[n-1]/a[n]*x^(n-1)+ … + a[1]/a[n]*x+a[0]/a[n]の形へ
double r=1.0, r0, t; //初期値を仮定する
for (i=1; i<=N; i++) {
r0=pow(fabs(N*a[i]),1.0/i);
if (r<r0) r=r0;
}
for (i=1; i<=N; i++) { //半径rの円に等間隔に配置する
t=2.0*PI/N*(i-3.0/4.0);
x1[i-1]=-a[N-1]/N+r*cos(t); //アーバスの初期値
x2[i-1]=r*sin(t);
}
double e,f1,f2,w1,w2,p1,p2,a1,a2,norm;
do {
e=0.0;
for (i=0; i<N; i++) {
f1=1.0; //分子 f(zi)、実部=f1,虚部=f2
f2=0.0; //※a(0)=1
for (j=N-1; j>=0; j--) { //ホーナー法 f(z)=( … ((z+a[n-1])*z+a[n-2])*z+a[n-3])* … +a[1])*z+a[0]
w1=f1*x1[i]-f2*x2[i]; //f*z=(f1+i*f2)*(x1+i*x2)=(f1*x1-f2*x2)+i*(f1*x2+f2*x1)
w2=f2*x1[i]+f1*x2[i];
f1=w1+a[j]; //f=z+a[j]
f2=w2;
}
p1=1.0; //分母 Π[j=1,N,i≠j](zi-zj)、実部=p1,虚部=p2
p2=0.0;
for (j=0; j<N; j++) {
if (j!=i) {
w1=p1*(x1[i]-x1[j])-p2*(x2[i]-x2[j]); //p*(zi-zj)
w2=p1*(x2[i]-x2[j])+p2*(x1[i]-x1[j]);
p1=w1; //p=(zi-zj)
p2=w2;
}
}
t=p1*p1+p2*p2; //分子÷分母 (f1+i*f2)/(p1+i*p2)=(f1+i*f2)(p1-i*p2)/(p1*p1+p2*p2)
if (t==0.0) {
printf("０では割れません。\n");
return 1;
}
a1=(f1*p1+f2*p2)/t;
a2=(f2*p1-f1*p2)/t;
norm=sqrt(a1*a1+a2*a2);
if (e<norm) e=norm; //最大値
x1[i]-=a1; //k回目の近似根 zi[k+1]=zi[k]-f(zi[k])/Π[j=1,N,i≠j](zi[k]-zj[k])
x2[i]-=a2;
}
} while (e>1.0e-10); //収束するまで ※調整が必要である
for (i=0; i<N; i++) printf( "x%d 実部=%f 虚部=%f\n", i,x1[i],x2[i]);
return 0;
}
実行結果
x0 実部=1.000000 虚部=0.000000
x1 実部=2.000000 虚部=-0.000000
x2 実部=3.000000 虚部=-0.000000
x3 実部=5.000000 虚部=-0.000000
x4 実部=8.000000 虚部=-0.000000
x5 実部=7.000000 虚部=-0.000000
x6 実部=6.000000 虚部=0.000000
x7 実部=4.000000 虚部=0.000000
続行するには何かキーを押してください . . .
|
|
|
|
投稿者：山中和義
投稿日：2014年10月17日(金)15時13分50秒
|
|
|
|
> No.3526[元記事へ]
tibitaさんへのお返事です。
complex.h を使ったサンプルです。 CPad 2.31 + Borland C++ Complier 5.5 で確認しました。
//DKA（Durand-Kerner-Aberth）法によるn次代数方程式f(x)=0の解法
#include <stdio.h>
#include <math.h>
#include <complex.h>
#define PI 3.14159265358979323846
#define N 8 //n次
int main(void)
{
double a[N+1] = {40320,-109584,118124,-67284,22449,-4536,546,-36,1}; //係数
complex<double> x[N]; //n個の根（実部、虚部）
int i,j;
for (i=0; i<=N; i++) a[i]/=a[N]; //x^n+a[n-1]/a[n]*x^(n-1)+ … + a[1]/a[n]*x+a[0]/a[n]の形へ
double r=1.0, r0; //初期値を仮定する
for (i=1; i<=N; i++) {
r0=pow(fabs(N*a[i]),1.0/i);
if (r<r0) r=r0;
}
for (i=0; i<N; i++) { //半径rの円に等間隔に配置する
x[i]=-a[N-1]/N+r*exp(2.0*PI/N*((i+1)-3.0/4.0)*(0,1)); //アーバスの初期値
}
double e,norm;
do {
e=0.0;
for (i=0; i<N; i++) {
complex<double> f,p,t;
f=1.0; //分子 f(zi) ※a(n)=1
for (j=N-1; j>=0; j--)
f=f*x[i]+a[j]; //ホーナー法 f(z)=( … ((z+a[n-1])*z+a[n-2])*z+a[n-3])* … +a[1])*z+a[0]
p=1.0; //分母 Π[j=1,N,i≠j](zi-zj)
for (j=0; j<N; j++)
if (j!=i) p*=x[i]-x[j]; //p*(zi-zj)
//分子÷分母 (f1+i*f2)/(p1+i*p2)=(f1+i*f2)(p1-i*p2)/(p1*p1+p2*p2)
if (abs(p)==0.0) {
printf("０では割れません。\n");
return 1;
}
t=f/p;
norm=abs(t);
if (e<norm) e=norm; //最大値
x[i]-=t; //k回目の近似根 zi[k+1]=zi[k]-f(zi[k])/Π[j=1,N,i≠j](zi[k]-zj[k])
}
} while (e>1.0e-10); //収束するまで ※調整が必要である
for (i=0; i<N; i++) printf("x%d 実部=%f 虚部=%f\n", i,real(x[i]),imag(x[i]));
return 0;
}
実行結果
D:\My Documents\C>test
x0 実部=3.000000 虚部=0.000000
x1 実部=4.000000 虚部=0.000000
x2 実部=7.000000 虚部=0.000000
x3 実部=2.000000 虚部=0.000000
x4 実部=1.000000 虚部=0.000000
x5 実部=6.000000 虚部=0.000000
x6 実部=5.000000 虚部=0.000000
x7 実部=8.000000 虚部=0.000000
-- Press any key to exit (Input "c" to continue) --
|
|
|
|
投稿者：山中和義
投稿日：2014年10月17日(金)16時33分7秒
|
|
|
|
> No.3529[元記事へ]
> //DKA（Durand-Kerner-Aberth）法によるn次代数方程式f(x)=0の解法
>
> #include <stdio.h>
> #include <math.h>
> #include <complex.h>
>
> #define PI 3.14159265358979323846
> #define N 8 //n次
>
> int main(void)
> {
> double a[N+1] = {40320,-109584,118124,-67284,22449,-4536,546,-36,1}; //係数
> complex<double> x[N]; //n個の根（実部、虚部）
> int i,j;
>
> for (i=0; i<=N; i++) a[i]/=a[N]; //x^n+a[n-1]/a[n]*x^(n-1)+ … + a[1]/a[n]*x+a[0]/a[n]の形へ
>
> double r=1.0, r0; //初期値を仮定する
> for (i=1; i<=N; i++) {
> r0=pow(fabs(N*a[i]),1.0/i);
> if (r<r0) r=r0;
> }
> for (i=0; i<N; i++) { //半径rの円に等間隔に配置する
> x[i]=-a[N-1]/N+r*exp(2.0*PI/N*((i+1)-3.0/4.0)*(0,1)); //アーバスの初期値
> }
>
> double e,norm;
> do {
> e=0.0;
>
訂正します。
int main(void)
{
double a[N+1] = {40320,-109584,118124,-67284,22449,-4536,546,-36,1}; //係数
complex<double> Zi(0.0,1.0); //虚数単位
complex<double> x[N]; //n個の根（実部、虚部）
int i,j;
for (i=0; i<=N; i++) a[i]/=a[N]; //x^n+a[n-1]/a[n]*x^(n-1)+ … + a[1]/a[n]*x+a[0]/a[n]の形へ
double r=1.0, r0; //初期値を仮定する
for (i=1; i<=N; i++) {
r0=pow(fabs(N*a[i]),1.0/i);
if (r<r0) r=r0;
}
for (i=0; i<N; i++) { //半径rの円に等間隔に配置する
x[i]=-a[N-1]/N+r*exp(2.0*PI/N*((i+1)-3.0/4.0)*Zi); //アーバスの初期値
}
double e,norm;
do {
e=0.0;
：
：
実行結果
D:\My Documents\C>test
x0 実部=1.000000 虚部=0.000000
x1 実部=2.000000 虚部=-0.000000
x2 実部=3.000000 虚部=-0.000000
x3 実部=5.000000 虚部=0.000000
x4 実部=8.000000 虚部=-0.000000
x5 実部=7.000000 虚部=0.000000
x6 実部=6.000000 虚部=-0.000000
x7 実部=4.000000 虚部=-0.000000
-- Press any key to exit (Input "c" to continue) --
|
|
|
戻る