求助

[wuximing]

已知：a R b

－674.843－1.54719×10^15x－2.13108×10^12ax＋7.10361×10^11aRx－2.80726×10^24x2－3.25852×10^24ax2＋7.41539×10^21b x2＋1.62862×10^24aR x2－4.3949×10^31x3－5.88422×10^33ax3＋1.63691×10^34bx3＋2.95501×10^33aRx3＋2.8072×10^38x4－4.64249×10^40ax4＋2.59434×10^41bx4＋4.6425×10^40aRx4＋4.41037×10^45x5＋3.42×10^43bx5＋6.12×10^42aRx5＋5.814×10^47x6=0
根据此方程求解X，X共有6个根，其中只有一个是正解。

请问有好点的算法吗?

[saint001]

一元6次方程可以求出6个复解
我用的是矩阵特征值的方法，求出所有的解，有程序

[saint001]

你给个例子，我算一下

[Linevan]

UP
好厉害的程序！

[saint001]

不是我写的
主要的三个函数是改编了了别人的

[HUNTON]

用牛顿迭代试试

[saint001]

牛顿迭代法只能求实根，并且对初值的依赖性很大，每次只能得到一个解，不知道什么时候、给多少个初值才能求出所有的解(实)
对于一般多项式求根问题，有精确的数值算法来求出所有的根，包括实根和复根，可以转化成矩阵的特征值求解，而特征值求解是有很多数值方法的。
n次方程an*x^n+a(n-1)*x^(n-1)+...+a0=0,(a(n-1)!=0)的所有复根等同于n*n矩阵
-a(m-1)/a(m) -a(m-2)/a(m) ... ... -a(1)/a(m) -a(0)/a(m)
1 0 ... ... 0 0
0 1 ... ... 0 0
............................................................
0 0 ... ... 1 0
的所有特征值
我有此方法的程序，试了一下你的方程，求的解都是实数，但不知道a,b,R怎样给才比较均衡

[saint001]

#include "math.h"
#include "stdio.h"

#define SIGN(a,b) ((b) > 0 ? fabs(a) : -fabs(a))
#define SWAP(g,h) {y=(g);(g)=(h);(h)=y;}
#define RADIX 2.0

void balanc(double **a, int n)
{
int last,j,i;
double s,r,g,f,c,sqrdx;
sqrdx=RADIX*RADIX;
last=0;
while (last == 0)
{
last=1;
for (i=1;i<=n;i++)
{
r=c=0.0;
for (j=1;j<=n;j++)
if (j != i)
{
c += fabs(a[j][i]);
r += fabs(a[i][j]);
}
if (c && r)
{
g=r/RADIX;
f=1.0;
s=c+r;
while (c<g)
{
f *= RADIX;
c *= sqrdx;
}
g=r*RADIX;
while (c>g)
{
f /= RADIX;
c /= sqrdx;
}
if ((c+r)/f < 0.95*s)
{
last=0;
g=1.0/f;
for (j=1;j<=n;j++)
a[i][j] *= g;
for (j=1;j<=n;j++)
a[j][i] *= f;
}
}
}
}
}

void elmhes(double **a,int n)
{
int m,j,i;
double y,x;

for (m=2;m<n;m++) {
x=0.0;
i=m;
for (j=m;j<=n;j++) {
if (fabs(a[j][m-1]) > fabs(x)) {
x=a[j][m-1];
i=j;
}
}
if (i != m) {
for (j=m-1;j<=n;j++) SWAP(a[i][j],a[m][j])
for (j=1;j<=n;j++) SWAP(a[j][i],a[j][m])
}
if (x) {
for (i=m+1;i<=n;i++) {
if (y=a[i][m-1]) {
y /= x;
a[i][m-1]=y;
for (j=m;j<=n;j++)
a[i][j] -= y*a[m][j];
for (j=1;j<=n;j++)
a[j][m] += y*a[j][i];
}
}
}
}
}

void hqr(double **a,int n,double *wr,double *wi)
{
int nn,m,l,k,j,its,i,mmin;
double z,y,x,w,v,u,t,s,r,q,p,anorm;

anorm=fabs(a[1][1]);
for (i=2;i<=n;i++)
for (j=(i-1);j<=n;j++)
anorm += fabs(a[i][j]);
nn=n;
t=0.0;
while (nn >= 1) {
its=0;
do {
for (l=nn;l>=2;l--) {
s=fabs(a[l-1][l-1])+fabs(a[l][l]);
if (s == 0.0) s=anorm;
if (fabs(a[l][l-1]) + s == s) break;
}
x=a[nn][nn];
if (l == nn) {
wr[nn]=x+t;
wi[nn--]=0.0;
} else {
y=a[nn-1][nn-1];
w=a[nn][nn-1]*a[nn-1][nn];
if (l == (nn-1)) {
p=0.5*(y-x);
q=p*p+w;
z=sqrt(fabs(q));
x += t;
if (q >= 0.0) {
z=p+SIGN(z,p);
wr[nn-1]=wr[nn]=x+z;
if (z) wr[nn]=x-w/z;
wi[nn-1]=wi[nn]=0.0;
} else {
wr[nn-1]=wr[nn]=x+p;
wi[nn-1]= -(wi[nn]=z);
}
nn -= 2;
} else {
if (its == 30) printf("Too many iterations in HQR");
if (its == 10 || its == 20) {
t += x;
for (i=1;i<=nn;i++) a[i][i] -= x;
s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
y=x=0.75*s;
w = -0.4375*s*s;
}
++its;
for (m=(nn-2);m>=l;m--) {
z=a[m][m];
r=x-z;
s=y-z;
p=(r*s-w)/a[m+1][m]+a[m][m+1];
q=a[m+1][m+1]-z-r-s;
r=a[m+2][m+1];
s=fabs(p)+fabs(q)+fabs(r);
p /= s;
q /= s;
r /= s;
if (m == l) break;
u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
if (u+v == v) break;
}
for (i=m+2;i<=nn;i++) {
a[i][i-2]=0.0;
if (i != (m+2)) a[i][i-3]=0.0;
}
for (k=m;k<=nn-1;k++) {
if (k != m) {
p=a[k][k-1];
q=a[k+1][k-1];
r=0.0;
if (k != (nn-1)) r=a[k+2][k-1];
if (x=fabs(p)+fabs(q)+fabs(r)) {
p /= x;
q /= x;
r /= x;
}
}
if (s=SIGN(sqrt(p*p+q*q+r*r),p)) {
if (k == m) {
if (l != m)
a[k][k-1] = -a[k][k-1];
} else
a[k][k-1] = -s*x;
p += s;
x=p/s;
y=q/s;
z=r/s;
q /= p;
r /= p;
for (j=k;j<=nn;j++) {
p=a[k][j]+q*a[k+1][j];
if (k != (nn-1)) {
p += r*a[k+2][j];
a[k+2][j] -= p*z;
}
a[k+1][j] -= p*y;
a[k][j] -= p*x;
}
mmin = nn<k+3 ? nn : k+3;
for (i=l;i<=mmin;i++) {
p=x*a[i][k]+y*a[i][k+1];
if (k != (nn-1)) {
p += z*a[i][k+2];
a[i][k+2] -= p*r;
}
a[i][k+1] -= p*q;
a[i][k] -= p;
}
}
}
}
}
} while (l < nn-1);
}
}

void eig(double **a,double *wr,double *wi,int n)
{
balanc(a,n);
elmhes(a,n);
hqr(a,n,wr,wi);
}

void solve(double *aa,double *wr,double *wi,int n)
//aa:1*(n+1),x:1:n
{
int i,j;
double **a;
a=new double*[n+1];
for(i=0;i<n+1;i++)
a[i]=new double[n+1];
for(i=0;i<n+1;i++)
for(j=0;j<n+1;j++)
a[i][j]=0;
for(i=1;i<n+1;i++)
a[1][i]=-aa[i]/aa[0];
for(i=1;i<n;i++)
a[i+1][i]=1;
eig(a,wr,wi,n);
}

int main1(int argc, char* argv[])
{
int i,j,n;
double **a,*wr,*wi;
FILE *fp=fopen("d:\\hqr.txt","r");
fscanf(fp,"%d",&n);
a=new double*[n+1];
wr=new double[n+1];
wi=new double[n+1];
for(i=0;i<n+1;i++)
a[i]=new double[n+1];
for(i=1;i<n+1;i++)
for(j=1;j<n+1;j++)
fscanf(fp,"%lf",a[i]+j);
fclose(fp);

printf("原矩阵\n");
for(i=1;i<n+1;i++)
{
for(j=1;j<n+1;j++)
printf("%f\t",a[i][j]);
printf("\n");
}
printf("\n");

eig(a,wr,wi,n);

printf("特征值\n");
for(i=1;i<n+1;i++)
printf("%f\t+i*%f\n",wr[i],wi[i]);
return 0;
}

void main()
{
int i;
const int n=6;
double A=1e-2,B=1e2,R=1e-2;
double aa[n+1],wr[n+1],wi[n+1];
//aa:方程系数，从高次到低次
aa[0]=5.814e47;
aa[1]=4.41037e45+3.42e43*B+6.12e42*A*R;
aa[2]=2.8072e38-4.64249e40*A+2.59434e41*B+4.6425e40*A*R;
aa[3]=-4.3949e31-5.88422e33*A+1.63691e34*B+2.95501e33*A*R;
aa[4]=-2.80726e24-3.25852e24*A+7.41539e21*B+1.62862e24*A*R;
aa[5]=1.54719e15-2.13108e12*A+7.10361e11*A*R;
aa[6]=-674.843;
solve(aa,wr,wi,n);
printf("方程的解\n");
for(i=1;i<n+1;i++)
printf("%18.16f\t+i*%18.16f\n",wr[i],wi[i]);
}